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Abstract 

The  Air  Breathing  Medium  Range  Air-to-Air  Missile  (ABMRAAM)  represents 
developmental  technology  which  incorporates  both  a  rocket  engine  and  a  RAMJET 
engine.  Such  a  missile  uses  proportional  navigation  guidance  plus  an  additional 
trajectory  loft  command.  This  thesis  examines  the  optimal  trajectory,  and  hence  the 
optimal  lofting  command,  for  an  ABMRAAM.  A  numerical  simulation  of  the  missile 
is  presented  and  the  necessary  conditions  for  the  optimal  trajectory  are  derived. 
From  these  conditions,  the  problem  can  be  numerically  solved  for  the  optimal  loft 
command. 
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Formulation  of  the  Optimal  Trajectory  for  an  Air-Breathing  Medium 

Range  Air-to-Air  Missile 


1.1  Background 


1.  Introduction 


The  United  States  Air  Force  (USAF)  is  investigating  the  feasibility  of  air-to- 
air  missiles  that  incorporate  both  a  rocket  engine  and  a  variable  throttle  RAMJET 
engine.  The  missile  with  both  a  rocket  engine  and  a  variable  throttle  RAMJET  has 
been  termed  a  Variable  Flow  Ducted  Rocket  (VFDR).  The  goal  is  to  improve  the 
performance  of  the  missile  for  the  same  preflight  weight  as  a  purely  rocket-propelled 


version.  Currently,  these  missiles  incorporate  proportional  navigation  guidance  plus 
a  loft  command,  gioft-  Generally,  the  loft  command  is  related  to  the  ratio  of  the 
current  line  of  sight  distance  to  the  original  line  of  sight  distance  between  the  missile 
and  the  target  and  an  altitude  command.  As  the  range  between  the  two  decreases, 
the  magnitude  of  the  command  decays  exponentially: 


9loft  —  Galtcmd 


(1.1) 


where 


gaitcmd  ~  altitude  command, 

Rlos  =  the  current  range  along  the  line  of  sight  vector, 
RloSo  =  the  initial  range  along  the  line  of  sight  vector,  and 
=  the  loft  decay  constant. 
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Due  to  the  short  duration  of  propulsion  and  near  constant  thrust,  commanding  a 
loft  according  to  Eq.  (1.1)  has  proven  to  be  effective  in  rocket-propelled  missiles 
as  a  method  of  increasing  the  potential  energy  of  the  vehicle.  Attempting  to  gain 
altitude,  thus  increasing  potential  energy,  may  or  may  not  be  the  most  effective 
guidance  implementation  for  a  variable  throttle  RAMJET.  Therefore,  the  optimal 
trajectory  of  this  missile  is  being  investigated  in  an  effort  to  determine  the  correct 
loft  command  for  the  VFDR  to  minimize  intercept  time. 

The  RAMJET  engine  compresses  the  incoming  air  without  the  mechanical 
compressor  associated  with  most  air-breathing  engines  (4).  The  pressure  of  the  in¬ 
coming  air  is  increased  by  passing  through  standing  shock  waves.  To  ensure  appro¬ 
priate  pressure  increase,  the  incoming  air  flow  must  remain  relatively  undisturbed. 
To  prevent  disturbances,  the  sideslip  angle  must  be  approximately  zero  for  the  VFDR 
due  to  the  configuration  of  the  RAMJET  inlet.  Hence,  the  missile  equipped  with 
a  RAMJET  must  fly  a  bank  to  turn  (BTT)  profile  to  maintain  zero  sideslip  angle. 
Therefore,  the  VFDR  utilizes  coordinated  turning  to  alter  its  cross  range  position, 
rather  than  the  skid  to  turn  (zero  roll  angle)  profile  of  the  typical  rocket  propelled 
missile. 

1.2  Previous  Work 

The  model  for  the  simulation  is  based  on  the  simulation  ENGAGE  (10).  EN¬ 
GAGE  was  created  for  use  with  personal  computers  and  simulates  a  one-on-one 
aircraft  pursuit  and  evasion.  Each  aircraft  is  capable  of  launching  one  missile  at  the 
other  aircraft.  The  algorithms  are  versatile  and  permit  extensive  variations  of  guid¬ 
ance  and  control  algorithms  for  both  the  aircraft  and  the  missiles  they  fire.  However, 
ENGAGE  is  not  useful  as  a  tool  to  optimize  a  given  performance  criterion;  whether 
the  criterion  is  minimum  time,  maximum  range,  etc.  Parameters  such  as  feedback 
gains,  loft  command  cutoff  values,  and  altitude  commands  can  be  varied  by  the  user 
in  an  attempt  to  achieve  improved  performance.  This  type  of  optimization  is  ad  hoc 
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at  best.  The  variable  throttle  RAMJET  complicates  performance  evaluation  even 
further  by  greatly  increasing  the  variation  of  the  missile’s  performance  based  upon 
the  trajectory  flown  (6:319).  Paris  et.  al.  evaluated  the  optimal  trajectory  for  an 
air-breathing  missile  problem  but  only  for  the  RAMJET  phase  of  boost.  ENGAGE 
simulates  the  full  flight  of  the  missile  but  is  not  designed  to  allow  the  application  of 
the  calculus  of  variations  to  determine  optimal  loft  command. 

1.3  Project  Scope 

The  scope  of  this  project  is  to  first  recreate  a  simulation  for  the  VFDR  using 
MATLAB™.  During  the  course  of  recreating  the  missile  portion  of  ENGAGE,  a 
few  errors  were  discovered  so  an  analysis  of  the  impact  is  presented.  Second,  the 
necessary  conditions  for  the  optimal  loft  control  are  developed.  The  minimum  time 
optimal  loft  control  for  the  three-dimensional  interception  problem  is  shown  to  be 
most  interesting  when  the  engagement  is  in  a  vertical  plane.  Therefore,  the  emphasis 
is  shifted  to  determining  the  optimal  control  assuming  missile  flight  is  restricted 
to  the  vertical  plane.  Third,  a  numerical  shooting  method  is  developed  to  solve 
the  two-point  boundary-value  problem  that  is  derived  as  a  result  of  optimizing  the 
loft  command.  The  discussion  so  far  has  frequently  referred  to  “optimal  control;” 
however,  optimal  performance  is  never  complete  until  a  measure  of  performance  is 
presented.  The  goal  of  this  analysis  is  to  find  the  optimal  control  for  loft  to  minimize 
the  final  time  of  rendezvous  with  the  target.  Hence,  the  performance  index  can  be 
stated  as 

minimize  J  =  tf  (1-2) 

where 


J  =  the  performance  index,  and 
tf  =  the  final  time. 
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1-4  Sequence  of  Presentation 

The  first  step  in  solving  this  problem  is  to  develop  a  simulation  for  the  VFDR. 
Chapters  two  and  three  describe  the  model  that  is  simulated  in  ENGAGE.  The 
description  includes  the  coordinate  systems,  the  atmospheric  model,  the  equations 
of  motion,  the  engines  and  the  guidance  and  control.  Chapter  four  analyzes  the 
simulations  and  discusses  the  errors  discovered.  Chapter  five  presents  the  derivation 
of  the  necessary  conditions  for  an  extremal  solution  for  the  free  final  time  problem. 
The  shooting  method  is  presented  as  a  numerical  technique  for  solving  the  two-point 
boundary- value  problem  formulated  from  the  necessary  conditions  for  an  extremal 
solution.  An  example  problem  and  the  numerical  results  accompany  the  discussion. 
The  formulation  for  the  optimal  trajectory  is  presented  for  the  three-dimensional 
intercept.  An  interesting  situation  occurs  when  flight  is  restricted  to  the  vertical 
plane.  Therefore,  the  two-dimensional  intercept  is  examined  and  the  formulation  for 
the  optimal  trajectory  two-point  boundary-value  problem  is  presented. 
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II.  Basic  Equations  for  the  Simulation 


2.1  Atmosphere 


The  atmosphere  is  modeled  using  the  U.S.  Standard  Atmosphere,  1962  (1). 
The  model  assumes  that  the  atmosphere  is  comprised  of  ideal  air  that  contains  no 
moisture  or  dust  and  can  be  characterized  by  the  ideal  gas  law.  An  ideal  gas  is  a 
gas  in  which  the  molecules  are  sufficiently  far  apart  so  that  intermolecular  forces 
are  negligible;  the  gas  acts  as  a  continuous  material  in  which  the  properties  are 
determined  by  statistical  average  of  the  particle  effects  (8:48).  Up  to  altitudes  of  80 
kilometers  (km),  air  is  assumed  to  be  homogeneously  mixed  with  a  relative  volume 
composition  that  leads  to  a  constant  molecular  weight,  MWaW-  Applying  the  ideal 
gas  law  to  air  yields 


pRuT 

MW~^ 


(2.1) 


where 


Ru  =  the  universal  gas  constant, 
P  =  the  total  pressure  of  air, 

T  =  the  total  temperature,  and 
p  =  density. 


However,  since  the  molecular  weight  of  air  is  assumed  constant. 


Ra 


Ru 


=  constant. 


Therefore,  the  ideal  gas  law  for  air  becomes, 

P  =  pRairT. 


(2.2) 
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To  80  km,  the  modeled  atmosphere  is  also  assumed  to  be  in  equilibrium  which  means 
that  the  pressure  is  related  to  the  geometric  altitude  by  the  differential  relation 

dP  =  —pgdho  (2-3) 


where 


ho  =  the  geometric  height  above  the  earth’s  surface,  and 
g  =  gravity  as  a  function  of  altitude. 


Equation  (2.3)  is  the  hydrostatic  equation.  The  customary  approach  in  standard 
atmosphere  calculations  is  to  eliminate  the  variation  of  gravity  with  altitude  from  the 
hydrostatic  equation.  To  effectively  eliminate  the  dependence,  geopotential  altitude 
h  is  introduced.  Before  the  geopotential  altitude  can  be  defined,  the  variation  of 
gravity  with  geometric  altitude  is  examined. 

Consider  two  small  elements  of  the  atmosphere,  where  one  is  just  at  the  surface 
of  the  earth  {ha  =  0)  and  the  other  is  at  some  geometric  altitude,  above  the 
surface  of  the  earth.  The  inverse  square  law  of  gravitation  expresses  the  attractive 
force  between  each  element  and  the  earth  as 


m^g  =  G 


m^g  =  G 


MeTTIi 

Re 

MEm2 

{Re  +  hoY 


(2.4) 

(2.5) 


where 


Re  =  the  radius  of  the  earth, 

G  =  the  gravitational  constant, 
mi  and  m2  =  the  mass  of  the  elements,  and 
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Me  —  the  mass  of  the  earth. 


Defining  g  as  the  acceleration  due  to  gravity  at  the  earth’s  surface,  it  is  noted  that 
g  =  g  when  ha  is  zero.  Also  noting  that  the  mass  of  each  element  cancels  from  the 
expressions  above,  Eqs.  (2.4)  and  (2.5)  become 


g  =  G- 


g  =  G 


{Re  +  hoY 


An  expression  for  5  as  a  function  of  geometric  altitude  is  obtained  by  forming  the 
ratio  of  the  above  equations, 

f = 


9=^9 


{Re  +  hoY 


The  hydrostatic  equation,  Eq.  (2.3),  can  be  expressed  either  as 


dP  =  —pgdho 


(2.10) 


or  noting  that  li  g  =  g,  then  dhc  is  equal  to  an  arbitrarily  small  change  in  geopo¬ 
tential  height,  dh,  and  the  hydrostatic  equation  can  be  expressed  as 


dP  =  —pgdh. 


(2.11) 


Equating  Eqs.  (2.10)  and  (2.11)  and  solving  for  dh  yields 


dh  =  -dhc- 
9 


(2.12) 
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Temperature  Profile  in  the  Standard  Atmosphere 


Figure  2.1  Standard  Atmosphere  Relation  Between  Temperature  and  Altitude 


Substituting  the  relation  obtained  in  Eq.  (2.8)  and  then  performing  the  integration 
an  expression  for  h  is 


dh  = 


h  = 


{Re  +  ha)^ 
Re 


dhc 


{Re  +  ho) 


he 


(2.13) 


Equation  (2.13)  is  the  geopotential  altitude,  an  altitude  that  assumes  gravity  is 
constant  and  equal  to  its  sea  level  value.  The  variation  of  temperature  with  altitude 
was  determined  experimentally  for  the  U.S.  Standard  Atmosphere,  1962  and  is  shown 
in  Fig.  2.1.  For  the  atmosphere  in  which  air-breathing  flight  takes  place,  there  are  two 
noticeable  traits.  There  is  a  temperature  gradient  region  and  an  isothermal  region. 
The  temperature  in  the  gradient  region  varies  linearly.  Therefore,  the  temperature 
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at  a  given  h  is  expressed  by 


r  =  T,  +  ~{h  ~  h) 


(2.14) 


where 


Tfc  =  the  reference  temperature  at  the  start  of  the  region  of  interest,  and 
hb  =  the  associated  reference  geopotential  altitude. 


Dividing  Eq.  (2.3)  by  Eq.  (2.1)and  performing  the  integration,  two  expressions  for 
density  are  formed.  The  first  expression  is  valid  for  the  region  that  temperature 
varies  linearly  with  h.  The  second  expression  is  applicable  in  the  isothermal  region. 


where. 


(2.15) 

f  n 

(2.16) 

Pb  =  the  reference  density  at  the  start  of  the  region  of  interest. 

Equations  (2.15)  and  (2.16)  are  two  expressions  for  p,  depending  upon  whether  the 
altitude  is  in  an  isothermal  region  or  a  temperature  gradient  region.  Table  2.1  gives 
the  corresponding  reference  values  for  h,  Tb,  pb,  and 

The  above  atmospheric  model  allows  the  determination  of  T,  p,  and  P  solely 
as  a  function  of  geopotential  altitude  up  to  80  km.  Figure  2.2  shows  the  variation 
of  density  and  pressure  as  a  function  of  geopotential  altitude  for  the  atmospheric 
model.  The  values  for  density  and  pressure  are  normalized  by  the  respective  sea- 
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Reference  Altitude, 
/i,  km 

Reference  Temperature, 

n,K 

Reference  Density, 
pb,  kg/m^ 

dT 

dh' 

K/km 

0 

288.15 

1.225 

-6.5 

11 

216.65 

0 

20 

216.65 

1 

32 

228.65 

1 

47 

270.65 

2.8 

52 

270.65 

0 

61 

252.65 

-2 

Table  2.1  Properties  of  the  Atmosphere  at  the  Isothermal  Gradient  Boundaries 


Density  and  Pressure,  Percent  Sea-Level  Value 


Figure  2.2  Density  and  Pressure  as  Percent  of  Sea  Level  Value 
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level  value.  Above  80  km,  air  begins  to  diffuse  and  the  molecular  weight  can  no 
longer  be  assumed  constant. 

From  these  quantities,  the  following  expressions  for  a,  the  speed  of  sound  in 
air,  and  Qo,  dynamic  pressure,  can  be  obtained: 


a  =  yj'TairRairT  (2.17) 

Qo  =  IfVl.  (2.18) 


'jair  =  the  the  specific  heat  ratio  of  air 
Vair  =  the  velocity  of  the  freestream  air. 


2.2  Coordinate  Systems 

This  section  describes  the  various  coordinate  systems  used  in  the  model  and 
the  necessary  rotation  sequences  used  to  define  the  axis  systems.  This  model  is 
developed  assuming  that  the  earth  is  a  fiat,  non-rotating  reference  frame.  This 
assumption  is  valid  given  the  relatively  short  duration  of  a  typical  flight  which  is  on 
the  order  of  one  minute. 

2.2.1  Inertial  Axis.  The  inertial  axis  is  formed  by  using  the  right  hand  rule 
with  xe  and  pe  forming  a  plane  parallel  to  the  surface  of  the  earth,  and  ze  normal 
to  the  surface  of  the  earth.  Figure  2.3  depicts  the  earth  fixed  coordinate  system. 

2.2.2  Velocity  Axis.  To  go  from  the  Earth-fixed,  or  inertial,  axis  system 
to  the  velocity  axis  system,  a  rotation  of  angle  ^  is  made  about  the  ze  axis.  This 
forms  an  intermediate  axis  xei^  Pei,  and  zei.  Figure  2.4  shows  this  rotation.  The 
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Z, 


Figure  2.3  Inertial  Axis  System 


Figure  2.4  Rotation  Sequence  From  Earth  to  Velocity  Axis 


transformation  from  the  inertial  axis  to  the  intermediate  axis  is  given  by 


'XEl' 

cos^  sin^  0 

'xe' 

VEX 

= 

—  sin  ^  cos  V’  0 

VE 

.zei  . 

0  0  1. 

-ZE. 

(2.19) 


Next,  a  rotation  of  angle  7  is  made  about  the  resulting  -yE\  axis.  This  is  chosen  so 
a  positive  7  represents  climbing  flight.  This  produces  unit  vectors  in  xv,  yv,  and  zy 


2-8 


Figure  2.5  Rotation  from  Velocity  to  Stability  Axis 


directions  with  the  total  velocity,  V,  always  pointed  in  the  xy  direction. 


'  Xy 

COS  7  0  sin  7’ 

xe\ 

yv 

= 

0  1  0 

Vei 

.zv . 

.  —  sin  7  0  cos  7 . 

.zei  . 

(2.20) 


Equations  (2.19)  and  (2.20)  can  be  combined  allowing  the  definition  of  Tev  to  be 
the  rotation  matrix  from  the  earth  to  the  velocity  coordinate  system. 


'  Xy ' 

cos  7  0  sin  7' 

cos^  sin^  O' 

Xe 

yv 

= 

0  1  0 

—  sin  ‘ip  cos  tp  0 

yE 

.Zy. 

_  —  sin  7  0  cos  7 . 

0  0  1. 

.ze  . 

cos  7  cos  ^  cos  7  sin  ^  sin  7' 

'xe' 

—  sin  ip  cos  xp  0 

yE 

.  —  sin  7  cos  xp  —  sin  7  sin  xp  cos  7  _ 

.  Ze  . 

Xe 


=  Tev 


VE 


i^E  A 


(2.21) 


(2.22) 


(2.23) 


2.2.3  Stability  Axis.  The  stability  axis  is  defined  by  a  roll  angle,  cr,  about 
the  -xv  axis.  This  results  in  coordinates  X5,  i/s,  and  zs-  Figure  2.5  shows  the 
rotation  sequence  to  the  stability  axis  system  and  is  given  by 
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Zs^i 

Figure  2.6  Rotation  Sequence  from  Stability  to  Body  Axis 


(2.24) 

(2.25) 


Equation  (2.24)  leads  to  the  definition  of  the  rotation  matrix  from  the  velocity  to 
the  stability  axis  system,  Tvs- 


2.2.4  Body  Axis.  The  body  axis  system  is  defined  relative  to  the  stability 
axis  system  by  two  rotations.  Figure  2.6  shows  the  rotation  sequence  from  the 
stability  axis  to  the  body  axis.  The  first  is  a  rotation  of  angle  /?  about  the  zs  axis. 
The  second  rotation  is  an  angle  a  about  the  resulting  -y  axis.  The  transformation 
from  the  stability  axis  to  the  body  axis  is 
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Figure  2.7  Body  Fixed  Axis  System 


'xfil  r  COSO  0  sin  a]  f  cos/3  sin^S  01 

ys  =  0  1  0  —  sm/3  cos^  0  ys  (2.26) 

.^B 1  L  —  sin  a  0  cos  aj[  0  0  lJL'2^5. 

cos  a  cos  /3  cos  a  sin  (3  sin  a  1  f  xs ' 

=  —  sin^  cos^  0  ys  (2.27) 

—  sin  a  cos  15  —  sin  a  sin  cos  a  J  [  . 

'xs' 

=  Tsb  ys  ■  (2.28) 

-ZS  . 

The  resulting  direction  vectors,  xg,  ys,  and  zb,  are  fixed  to  the  missile  with  xb  out 
the  nose,  yB  out  the  left  side  of  the  missile  and  zb  out  the  top;  see  Fig  2.7.  The 
transformation  from  the  stability  axis  system  to  the  body  axis  system  is  accomplished 
by  Tsb,  which  is  defined  by  Eq.  (2.28). 

2.3  Force  Equations 

Figure  2.8  is  a  free  body  diagram  of  the  missile  in  flight.  The  missile  is  assumed 
to  be  a  point  mass.  The  sum  of  forces  produces  the  net  acceleration  of  the  vehicle. 
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Axial 


X 


s.v 


Velocity 


Y  Weight 

Figure  2.8  Free  Body  Diagram  of  the  Missile  in  Flight 
The  acceleration  terms  are  derived  followed  by  a  development  of  the  force  terms. 


where 


=  may 


m 


^(Vv)  +  <^v/E  X  Vv 


(2.29) 


Fv  =  the  total  force  on  the  missile, 
ay  =  the  total  acceleration  of  the  missile, 

Vy  =  the  velocity  vector  in  the  velocity  frame, 

U)v/E  =  the  angular  velocity  vector  of  the  velocity  frame 
relative  to  the  inertial  frame. 
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The  forces  acting  on  the  missile  in  flight  are  in  various  reference  frames.  Weight,  W, 
is  acting  in  the  inertial  frame  and  the  axial,  Fa,  normal,  lateral,  Fy,  and  thrust, 
Fti  forces  are  in  the  body  frame. 

By  definition,  the  velocity  axis  system  is  aligned  with  the  velocity  vector. 


Vv 


'V 

0 

0 


(2.30) 


The  angular  velocity  of  the  velocity  axis  system  relative  to  the  inertial  frame  is 
known  through  the  rotations  required  to  go  from  the  inertial  to  velocity  coordinate 
systems. 

<^v/E  =  4’^e  —  iyE\  (2.31) 

However,  ipZE  and  ^yEi  must  be  expressed  in  the  velocity  frame,  hence 


and 


i’v  =  TeV 


0 

0 


cos  7  cos  xl>  cos  7  cos  ^  sin  7 ' 

’O' 

—  sin  t/}  cos  ^  0 

0 

,  —  sin  7  cos  V’  —  sin  7  sin  V’  cos  7  _ 

^  sin  7 
0 

^cos7 


cos  7  0  sin  7' 

’O' 

iv  = 

0  1  0 

7 

.  —  sin  7  0  cos  7 . 

.0. 

(2.32) 
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=  7 


Therefore,  can  be  expressed  entirely  in  the  velocity  axis  system  as 


'^vjE  =  i’  sin  'yxv  —  jyv  +  ^  cos  'yzy. 


(2.33) 


Knowing  ojv/e  and  Vvt  the  cross  product  of  Eq.  (2.29)  becomes 


XV  yv 


(jJv/e'X^v  =  '0sin7  —7  ipcos'f 


E^cos7 

V-y 


(2.34) 


Also, 


-(Vv)=  0 


=  Vxv- 


(2.35) 


Combining  Eqs.  (2.34)  and  (2.35)  the  total  acceleration  vector,  av  becomes 


av  =  ^(^v)  +  Wy/EX  Vv 


=  E^cos7  . 


(2.36) 
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For  completion  of  Eq.  (2.29),  the  net  forces  in  the  velocity  axis  frame  must  be 
determined.  The  forces  acting  in  the  body  frame  are  the  thrust,  the  axial  force,  the 
normal  force  and  the  side  force.  The  thrust  force  is  aligned  with  the  xb  axis.  The 
axial  force  is  also  aligned  with  the  xb  axis  but  acts  in  the  opposite  direction  of  the 
thrust  force.  The  normal  force  acts  normal  to  the  xb  —  Vb  plane  in  the  zb  direction. 
The  lateral  force  acts  in  the  yB  direction.  The  force  vector  is  expressed  as 


Fb 


Ft  —  Fa 
F 

^  y 
F 

L  •*  n  J 


(2.37) 


where  Fb  indicates  the  force  vector  is  only  thrust  and  aerodynamic  forces,  and 
gravity  is  not  yet  included.  As  a  first  step  in  transforming  Fb,  the  force  vector  is 
rotated  to  the  stability  axis. 


Fs  =  TbsFb 
==  Tsb^Fb 


cos  a  cos  0 

—  sin  0 

—  sin  oc  cos  0 ' 

'Ft  —  Fa' 

cos  a  sin  0 

cos  0 

—  sin  a  sin  0 

Fy 

sin  a 

0 

cos  a 

.  Fn  . 

Ft  cos  a  cos  ^  —  Fa  cos  a  cos  ^  —  Fy  sin  ^  —  Fn  sin  a  cos  {3 
Ft  cos  a  sin  /3  —  Fa  cos  a  sin  ^  +  Fy  cos  ^  —  Fn  sin  a  sin  /7 
Ft  sin  a  —  Fa  sin  a  +  Fn  cos  a 


(2.38) 


Equation  (2.38)  represents  the  thrust  and  aerodynamic  forces  acting  on  the  missile 
relative  to  the  stability  axis  system.  Drag,  D,  side,  S  and  lift,  T,  are  defined  in 
the  missile  stability  axis  and  are  pictured  in  Fig  2.9.  Lift  is  the  total  aerodynamic 
force  acting  on  the  missile  in  the  stability  axis  in  the  zs  direction.  Drag  is  the  total 
aerodynamic  force  acting  on  the  missile  in  the  stability  axis  in  the  -xs  direction.  The 
side  force  is  the  total  aerodynamic  force  acting  on  the  missile  in  the  stability  axis  in 


2-15 


the  ys  direction. 


Fa  COS  a  cos  ^  +  Fy  sin  0  +  Fn  sin  a  cos  ^ 

(2.39) 

—Fa  cos  a  sin  ^  +  Fy  cos  ^  —  Fn  sin  a  sin  ^ 

(2.40) 

—Fa  sin  a  -f  cos  a 

(2.41) 

The  above  definitions  reduce  Eq.  (2.38)  to 


Fs 


Ft  cos  a  cos  P  —  D 
Ft  cos  a  sin  ^  +  S 
Ft  sin  a L 


(2.42) 


The  final  rotation  sequence  from  the  stability  axis  system  to  the  velocity  axis  system 
is  accomplished  by  the  transpose  of  the  rotation  matrix  in  Eq.  (2.25).  The  result  of 
the  rotation  is  the  thrust  and  the  aerodynamic  forces  are  expressed  in  the  velocity 
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axis  system. 


■£>  _  rp  r-ri 

rv  —  J-vs 


1 

0 

0  ■ 

’  Ft  cos  a  cos  (3  —  D' 

0 

cos  a 

sincr 

Ft  cos  a  sin  ^  -1-  S' 

0 

—  sin  a 

cos  (7 

Ft  sin  a  -f  Z 

F-r  cos  a  cos  ^  —  D 


(Ft  sin  a  +  Z)  sin  <7  +  {Ft  cos  a  sin  ^  +  S)  cos  a 
{Ft  sin  a  +  L)  cos  a  —  {Ft  cos  a  sin  +  S')  sin  cr 


(2.43) 


Fv  is  the  thrust  and  aerodynamic  forces  acting  of  the  missile  in  the  velocity  reference 
frame.  Gravity  effects  are  now  included  to  produce  the  total  force  vector  in  the 
velocity  coordinate  system,  Fy.  Gravity  acts  in  the  -ze  direction  of  the  inertial 
frame.  Noting  the  mass  of  the  missile  is  m  and  using  the  relation 


We  =  rn 


■  0  ■ 
0 


L-S^J 


(2.44) 


yields 


We  = 


0  ■ 
0 

-w 


Resolving  the  weight  vector  into  the  the  velocity  axis  components  yields, 


Wv  =  Tev 


0  ■ 
0 

-w 


cos  7  cos  ^  cos  7  cos  ^  sin  7 ' 

■  0  ■ 

—  sin  ^  cos  Ip  0 

0 

_  —  sin  7  cos  Ip  —  sin  7  sin  ip  cos  7  _ 

-w_ 
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~W  sin  7 

0 

—W  cos  7 


(2.45) 


The  thrust  and  aerodynamic  forces,  Eq.  (2.43),  are  combined  with  the  weight, 
Eq.  (2.45)  to  denote  the  total  forces  acting  on  the  missile  in  the  velocity  frame. 


Fv  =  Fv  4" 


Ft  cos  a  cos  ^  —  D 

’  —  W  sin7  ’ 

{Ft  sin  a  +  T)  sin  cr  +  {Ft  cos  a  sin  +  S)  cos  a 

+ 

0 

{Ft  sin  a  L)  cos  cr  —  {Ft  cos  a  sin  /3  +  S)  sin  a . 

_—W  C0S7  J 

Ft  cos  a  cos  ^  —  D  —  W  sin  7 

{Ft  sin  a  +  L)  sin  <T  +  {Ft  cos  a  sin  ^  +  S')  cos  a 
{Ft  sin  a  +  L)  cos  a  —  {Ft  cos  asm/3  +  S)  sin  a  —  W  cos  7 


(2.46) 


The  expressions  for  the  acceleration  of  the  missile  in  the  velocity  axis  and 
the  total  forces  acting  on  the  missile  in  the  velocity  axis  can  now  be  combined  in 
Eq.  (2.29)  to  complete  the  expressions  for  the  equations  of  motion  of  the  missile. 
Recall  Eq.  (2.29); 


^Fv  =  mav 


Substituting  the  expression  obtained  from  Eq.  (2.46)  for  Fv  and  the  expression 
obtained  from  Eq.  (2.36)  for  ay  yields 


■  y  ■ 

Ft  cos  a  cos ^  —  D  —  W  sm^ 

m 

V  '0  cos  7 

= 

{Ft  sin  a  +  L)  sin  <7  +  {Ft  cos  a  sin  ^  +  5)  cos  cr 

^7 

.  {Ft  sin  a  +  T)  cos  a  —  {Ft  cos  q  sin  +  S')  sin  cr  —  VE  cos  7 . 

(2.47) 
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Dividing  both  sides  of  the  equation  by  the  mass  of  the  missile,  the  above  equation 
becomes 


V'lpcos'y  =  — 


Substituting 


Ft  cos  a  cos  ^  —  D  —  W  sin-y 

{Ft  sin  a  +  L)  sin  <t  +  {Ft  cos  a  sin  +  S)  cos  a 

{Ft  sin  a  +  L)  cos  a  —  {Ft  cos  a  sin  ^  +  5)  sin  cr  —  VF  cos  7 

(2.48) 


—  =  A 

m  W 


(2.49) 


into  Eq.  (2.48),  the  expression  becomes 

V  1  r  cos  a  cos  /3  -  -  sin  7) 

V^cos7  =  5  [(^  sina  + sin<T  +  cosasin^  + COSO- 

.  Vi  J  [(^sina  +  ^)coso-- (^cosasin/?  +  ^)sino--cos7 

Rearranging  the  terms,  expressions  for  y,  ip  and  7  are 

'VI  r  5^  cos  a  cos  -sin  7) 

^  [(^sina  +  -^)sino-+ (^cosasin^  +  -|:)cosa]  .  I 

.  7  J  L  [{^  sin  a  +  cos  cr  -  cos  a  sin  ^  sin  o-  -  cos  7  . 

Noting  5  =  0  and  ^  =  0  for  this  missile,  the  equations  of  motion  reduce  to 


(2.51) 


^1  r  S' (^cosa  - -  sin7) 

^  =  v^(l^sina  +  i)sin(r 

.  tJ  L  V  [(^  sina  + COSO- -  COS7 


(2.52) 


isi  VEi  and  ze  are  related  to  the  velocity  vector  in  the  velocity  axis  system 
through  the  transpose  of  Eq.  (2.23) 


ifE  =  Tev  0 
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'  V  cos  7  cos  ^ ' 

=  1^008  78111^  .  (2.53) 

y  sin  7 

Equation  (2.53)  gives  the  translational  equations  needed  to  determine  the  po¬ 
sition  of  the  missile  in  the  inertial  coordinate  system  and  Eq.  (2.52)  represents  the 
equations  of  motion  in  the  velocity  frame  for  the  missile  in  flight.  Specific  expressions 
for  Ft,  L  and  D  are  developed  in  the  following  sections. 

2.4  Aerodynamic  Forces 

The  aerodynamic  forces  for  this  missile  are  in  the  form  of  dimensionless  aero¬ 
dynamic  coefficients.  This  data  is  in  tabular  form  and  was  generated  through  wind 
tunnel  testing.  There  is  no  sideslip  angle,  hence  there  is  no  side  force.  The  axial 
force  is  modeled  as  a  function  of  a  dimensionless  axial  force  coefficient,  Cq,  based 
on  Mach  number,  M,  and  a.  A  dimensionless  compensation  to  the  axial  force  coef¬ 
ficient,  ACa,  as  a  function  of  altitude  and  Mach  number  corrects  for  the  variations 
in  aerodynamic  effects  as  density  and  Mach  number  vary.  The  axial  force.  Fa,  is 
expressed  as 


Fa  —  {Ca  -\-  ACa) Aref  Qo 

(2.54) 

Ca  =  f{a,M) 

(2.55) 

ACa  =  fih,M) 

(2.56) 

where 

Aref  =  the  aerodynamic  reference  area  of  the  missile. 

Two  sets  of  data  are  required  for  Ca.  The  RAMJET  requires  air  flow  into  the  com¬ 
bustion  chamber,  which  requires  ducts  to  open.  These  ducts  are  not  deployed  during 
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the  rocket  boost  phase,  hence  the  two  different  axial  force  profiles,  one  corresponding 
to  the  axial  coefficient  during  the  rocket  boost  phase,  and  the  other  corresponding 
to  the  axial  coefficient  during  the  RAMJET  boost  phase. 

The  normal  force,  is  modeled  as  a  function  of  a  dimensionless  normal  force 
coefficient,  C„,  based  on  M  and  a. 

Fn  =  CnArefQo  (2.57) 

C'„  =  /(a,M)  (2.58) 


2.5  Engine 

The  engine  needs  to  be  described  during  the  three  phases  of  its  flight.  Overall, 
there  are  two  propulsion  phases  and  one  coast  phase.  The  rocket  boost  phase  is  the 
first  phase.  The  rocket  booster  burns  at  a  fixed  rate  until  the  propellant  is  depleted. 
The  next  propulsion  phase  is  the  air-breathing  boost  phase.  This  phase  is  terminated 
when  the  oxidizer  is  depleted.  Unlike  the  constant  fuel  flow  rate  in  the  rocket  boost 
phase,  the  fuel  flow  rate  is  variable  in  the  air-breathing  boost  phase.  Separating  the 
rocket  boost  phase  and  the  air-breathing  boost  phase  is  a  brief  (0.2  sec)  coast  phase 
where  there  is  no  thrust  developed  and  no  fuel  depletion.  This  phase  begins  when 
the  propellant  for  the  rocket  engine  is  depleted  and  the  transition  time  permits  the 
opening  of  the  RAMJET  inlet.  Finally,  the  missile  returns  to  the  coast  phase  if  the 
oxidizer  is  depleted  prior  to  target  intercept. 

2.5.1  Rocket  Boost  Phase.  The  thrust,  Fj,  of  a  rocket  engine  pictured  in 
Fig.  2.10  is  given  by 

Ft  =  Tsi  +  AexitiPsi  -  p),  (2.59) 

where 


Tsi  =  the  thrust  at  sea  level. 
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Thrust,  Fj. 

- 


Figure  2.10  Rocket  Engine  Control  Volume 


Aexit  =  the  exit  area  of  the  nozzle, 

Psi  =  the  atmospheric  pressure  at  sea  level. 

This  equation  is  useful  since  T^i  is  known  from  test  data  and  Aexit  is  known  from  the 
geometry  of  the  missile.  The  pressure,  is  known  as  a  function  of  altitude  from 
the  atmospheric  model  and  is  given  by  Eq.  (2.1).  The  weight  of  the  missile,  Wm, 
decreases  at  the  rate  the  fuel  is  depleted,  Wf] 

w,  =  ^  (2.60) 

-•sPsI 

where 


Isp^^  =  the  sea  level  specific  impulse 

and  is  known  from  test  data.  Tgi  and  are  modeled  as  constant  values  which 
implies  that  wj  is  constant  for  the  rocket  boost  phase  of  flight. 

2.5.2  Air-Breathing  Boost  Phase.  A  general  RAMJET  engine  is  pictured 
in  Fig.  2.11  (11:96).  Air  at  station  0  enters  the  engine  at  free-stream  velocity,  Vo,  and 
pressure,  po,  at  a  rate  of  Wa.  The  capture  area  is  Aq.  The  velocity  of  the  entering  air 
is  reduced  and  the  static  pressure  is  increased  by  the  supersonic  diffuser  at  station 
1.  The  subsonic  diffuser,  station  2,  then  compresses  the  air  further.  The  air  flows 
into  the  combustor  at  station  3,  which  houses  the  burners.  The  air  is  heated  by  the 
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Inlet 

Figure  2.11  Features  of  the  General  RAMJET  Engine 

continuous  combustion  of  fuel.  The  heated  products  of  combustion  are  expanded  in 
the  nozzle,  station  4,  and  are  ejected  at  station  6  at  a  speed  greater  than  the  incoming 
air.  The  increase  in  momentum  of  the  gas  results  in  a  thrust  in  the  direction  of  flight. 

Since  this  engine  requires  flowing  air  in  order  to  build  pressure  for  the  combus¬ 
tion  process,  the  engine  is  not  able  to  generate  thrust  at  zero  flight  speed.  Hence,  the 
RAMJET  must  be  propelled  to  a  minimum  velocity  to  operate.  The  ABMRAAM 
accomplishes  this  through  the  rocket  booster  previously  described.  Once  in  flight, 
the  thrust,  Fy,  generated  by  the  RAMJET  is  given  by 

Ft  =  +  {Ve.it  -  Po)A,,iu  (2.61) 

9  9 

where 

Vexit  —  the  velocity  of  the  gases  exiting  the  nozzle, 

Pexit  =  the  pressure  of  the  gases  exiting  the  nozzle,  and 
Aexit  —  the  exit  area  of  the  nozzle. 
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A  throttle  is  used  to  control  the  fuel  flow  rate,  wj,  entering  the  combustion 
chamber  which  regulates  the  thrust.  The  throttle  logic  will  be  discussed  further 
in  the  next  chapter.  Once  Wj  is  determined,  Wa,  Vexit-,  and  Pexit  in  Eq.  (2.61)  are 
functions  of  M  and  a.  Again,  the  data  for  this  model  is  in  tabular  form. 

2.5.3  Coast  Phase.  As  previously  mentioned,  the  coast  phase  can  occur 
twice  during  a  missile  flight.  The  first  time  it  occurs  is  between  the  rocket  and  the 
RAMJET  boost  phases.  The  duration  is  for  0.2  seconds  and  allows  for  the  transition 
of  power  modes.  The  next  time  the  coast  phase  can  occur  is  if  intercept  has  not 
happened  prior  to  the  expenditure  of  the  fuel  to  sustain  the  RAMJET  propulsion. 
In  either  situation,  the  drag  induced  by  the  aerodynamic  effects  over  the  inlet,  Di  of 
the  RAMJET  is  given  by 

Di  =  Cd.QoAo  (2.62) 

where 


Cdi  =  the  coefficient  of  induced  drag. 


Cdi  is  a  function  of  M  and  a  and  is  in  tabular  form  for  the  missile  model.  This  table 
was  generated  by  wind  tunnel  testing. 
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III.  Guidance  and  Control 


3.1  Guidance 


The  VFDR  uses  a  mixed  guidance  strategy  to  guide  it  to  the  target.  The 
lateral  guidance  is  a  proportional  navigation  command.  The  horizontal  guidance  is 
a  combination  of  a  Mach  command  and  a  gravity  command.  The  vertical  guidance  is 
a  combination  of  a  loft  command,  a  proportional  navigation  command  and  a  gravity 
command.  This  section  will  describe  the  guidance  of  the  VFDR  and  derive  the 
line  of  sight  parameters  for  proportional  navigation.  Commands  not  generated  by 
proportional  navigation  will  then  be  discussed.  Proportional  navigation  guidance 
results  in  an  acceleration  command  in  the  velocity  axis  system  proportional  to  the 
angular  rate  of  change  of  the  line-of  sight  vector,  from  the  missile  to  the 

target.  N  is  called  the  proportional  navigation  constant.  Therefore,  the  commands 
generated  to  guide  the  missile  are  proportional  to  the  angular  velocity  of  the  missile 
relative  to  the  target  by  a  factor  of  N.  Generally,  N  ranges  from  2  to  6  which 
implies  that  the  missile  develops  a  lead  angle  on  the  target  (9).  If  N  =  I,  then  the 
resulting  missile  accelerations  will  alter  the  relative  velocity  between  the  missile  and 
the  target  and  drive  the  line  of  sight  rate  to  zero  and  if  iV  <  1,  the  missile  will  lag 
behind  the  target  (3:262).  Differentiating  the  velocity  vector  with  respect  to  time, 
Eq.  (3.1),  yields  the  desired  acceleration  vector,  a™,  of  the  missile  in  terms  of  the 
angular  velocity  of  line  of  sight  vector. 


3-in 


dt 


+  Nu^los  X  V 


=  Vxv  +  Nu>los  X  V. 


(3.1) 
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LOS 


Figure  3.1  Typical  Missile-Target  Intercept  Scenario 


Taking  the  cross  product  of  the  final  term  in  Eq.  (3.1)  yields 


Xv 

yv 

Zv 

0 

Nulos  X  V  = 

NUy 

NuJz 

= 

VmNiVz 

Vm 

0 

0 

.  Vm  . 

=  Vm^v  +  VmNUzyv  —  VmNuJyZv  (3.2) 


Normalizing  the  acceleration  vector,  am,  with  respect  to  gravity  yields 

(3 

Specific  expressions  for  Uy  and  in  Eq.  (3.3)  are  developed  by  examining  the  line  of 
sight  parameters  between  the  missile  and  the  target.  Note  that  does  not  appear  in 
Eq.  (3.2).  Figure  3.1  shows  a  general  three-dimensional  missile  target  game.  Before 
expressions  can  be  obtained  for  Uy  and  in  Eq.  (3.2),  some  initial  definitions  must 
be  made  for  the  components  of  distance  from  the  target  to  the  missile  and  the 


9 

Vm  .  ,  VmNuJz  .  V^NUy  , 

=  — XV  -f- - yv - 

9  9  9 


3-2 


components  of  relative  velocity  from  the  missile  to  the  target. 


^LOS  —  (^■^) 

Vlos  ^VT  —  ym  (3.5) 

zlos  ^hx  —  hm  (3.6) 

^^LOS  ~  COS  l/^m  COS  'Yjjj  (3.7) 

KiOS  =yT-Vm  sin  i>rn  COS  Jrn  (3.8) 

^Los  ~  'I'm  sin  (3.9) 


where  the  subscript  T  denotes  the  target  and  m  denotes  the  missile.  From  the 
definitions  in  Eqs  (3.4)  through  (3.9),  the  range  along  the  line  of  sight,  RloSi  can 
be  expressed  as 

1|72los||  =  (a^ios  +  y\os  +  ^ios)*  •  (3.10) 

The  closing  velocity,  Vdose  can  be  expressed  as 


^close 


Rlos  •  Vlos 

ll^iosll 

[  xlos  yios  Zlos  ]  •  [  V^los  ^ilos  ^L05  . 
Il^iosll 


(3.11) 


Dividing  the  line  of  sight  range  by  the  closing  velocity,  an  estimation  for  the  time 
remaining  until  intercept,  tgo  is  expressed  as 


tgo  — 


\\R^ 


LOS\ 


V. 


close 


(3.12) 


Referring  to  Fig  3.1,  the  rate  of  change  of  the  position  vector  from  the  missile  to 
the  target  can  be  described  by  a  cross  range  angular  velocity,  if^LOS,  and  an  angular 
rate  of  change  from  the  horizontal,  jias-  The  total  angular  velocity  of  the  line  of 
sight  vector  is 

LOS  =  IpLOSZE  —  iLOSysi-  (3.13) 
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^Ylos 


Figure  3.2  Line  of  Sight  Geometry  for  the  xlos  —  Vlos  Plane 

Specific  expressions  for  ipios  and  'Jlos  in  terms  of  the  definitions  for  the  line  of  sight 
parameters  are  required  to  generate  the  necessary  guidance  commands  in  Eq.  (3.3) 
to  intercept  the  target.  Examining  the  xios  —  Vlos  plane,  ipLOS  is  shown  in  Fig  3.2 
along  with  the  tangential  and  normal  components  of  the  line  of  sight  velocity  vector 
in  the  plane.  and  Vy^os  li^®  components  in  the  xlos  and 

Vlos  directions,  respectively.  Resolving  the  velocity  vectors  into  their  xei  and  yEi 
components  yields 

'V^^osE\'\  r  cos x/^LOs  sin V’LOs  01  \ Volos' 

^vlosEi  =  -  sin  xpLOS  cos  V’los  0  VyLos 

.14i05El-  -  0  0  1.  .V^LOS. 

COS  V’L05  +  sin  ^iOS  ' 

=  -14iosSinV’LOS  + Vj,^05  cos  V’ios  •  (3.14) 

V 

L  '^ULOS 

The  line  of  sight  velocity  component  in  the  yEi  direction  is  the  tangential  velocity 
of  point  A  that  is  being  swept  by  the  radius  {xIqs  +  ylos)^  with  angular  velocity 


3-4 


i^LOS-  Equating  the  two  expressions  for  the  velocity  of  point  A  yields 


(4os  +  ylos) '  i^LOS  =  cos  tl^LOS  -  Ki05  sin  ^ws-  (3.15) 


Therefore,  ij^LOS  can  be  expressed  as 


i’Los  = 


^VLos  cos  V>L05  -  Vxi^os  sin  i’LOS 

(Ax)S  +  vlos)^ 


(3.16) 


Noting  from  the  geometry  in  Fig  3.2  that 


,  XLOS 

cos  ip  LOS  = - r 

{^\os  +  VlosY 

.  ,  VLOS 

Sin  ipLOS  = - r 

(^io5  +  VlosY 


the  expression  for  ij^ios  becomes 


1  Yilos^I'OS  YciosVLOS 

VLOS  =  - 2 - ; - 2 - • 

^LOS  +  Vlos 


(3.17) 


The  rotation  due  to  the  7los  is  shown  in  Fig  3.3. 


^==LOsV 


''VLOSV  = 


^^LOSV  . 


■  COS7  0  sin7l  r  cos  V’los  +  l^ios  sin ' 

0  1  0  sin  ^iosT  14^05  cos  0LOS 

,-sin7  0  C0S7J  [ 

(Ki05  cos  ^lo5  +  sin  t/’Los)  cos  7  +  14^05  sin  7 

-VxLos  sin  V’L05  +  14i05  cos  if; LOS  (3.18) 

.  -  (Kto5  cos  IpLOS  +  sin  ipLOs)  sin  7  +  cos  7 . 


As  before,  the  tangential  components  of  velocity  are  equated 


RlosIlos  =  14ios  cos  7  -  (14x,os  IpLOS  +  V^j/ios  sin  ipLos)  sin  7.  (3.19) 
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Recall  the  angular  line  of  sight  velocity  vector 


<^LOS  —  i’LOSZE  —  iLOSyEl 

Rotating  i/^los^e  to  the  intermediate  coordinate  system  yields 


Xei 

cos^  sinV’  O' 

■  0  ■ 

Vei 

= 

—  sin  V’  cos  V’  0 

0 

.ze\  . 

0  0  1. 

.  4’los  . 

■  0  ■ 
0 

.  '4’los  . 


The  line  of  sight  angular  velocity  rate  becomes 

‘^losEI  =  fpLOSZEl  —  ihOSyEl- 


(3.22) 


(3.23) 


(3.24) 


Rotating  to  the  velocity  axis  system  produces 


COS  7  0  sin  7' 

■  0  ■ 

^LOS^  ~~ 

0  1  0 

— 7LOS 

.  —  sin  7  0  cos  7 . 

.  4los  . 

’  i’LOS  sin  7  ■ 
—iios 

.V’LOS  cos  7. 


(3.25) 


The  angular  rate  of  change  in  the  velocity  axis  of  the  line  of  sight  vector  is  given  by 
Eq.  (3.25).  Therefore, 


Wy  =  — 7los  (3.26) 

=  ikhos  cos  j  (3.27) 
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are  the  expressions  required  to  generate  a  proportional  navigation  command.  Recall 
Eq.  (3.3)  for  the  proportional  navigation  command  vector 


Vm  .  ,  V^NlO^  ^  VmNiJy  ^ 

rim  =  - d - yv - ^v- 

9  9  9 


(3.28) 


Substituting  Eqs.  (3.26)  and  (3.27)  yields 


Vm  .  ,  VmNlf^LOS  cos  ^  ^  VmN'^LOS. 

— XV  + - yv  + - 

9  9  9 


(3.29) 


which  is  the  dimensionless  command  vector.  The  proportional  navigation  command 


vector  used  is 


dcmdP/N  - 


VmNll>i,nf;  cos 'i 


VmAT'Vf.O.c; 


(3.30) 


Proportional  navigation  is  not  the  only  command  used  to  guide  the  missile. 
The  other  commands  employed  by  the  VFDR  will  be  examined.  Only  proportional 
navigation  is  used  in  the  lateral,  direction.  Therefore,  the  total  gcmdy  is  given  by 


9cmdy  —  9cmdP/ny  — 


Vm  N  Ip  LOS  COS  J 


(3.31) 


where 


9cmdy  =  the  total  command  in  the  yv  direction. 


(3.32) 


The  command  in  the  horizontal  direction  is  a  combination  of  a  command  to  control 
the  Mach  number  and  a  command  to  counter  the  effects  of  gravity.  The  Mach 
command,  gmachx  is  generated  first  comparing  a  desired  Mach  number  to  the  actual 
Mach  number  of  the  missile; 


Mcmd  -M  =  AM 


(3.33) 
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where 


Mcmd  =  desired  Mach  number 
M  =  actual  Mach  number 
AM  =  delta  Mach  number,  desired  -  actual. 


The  VFDR  uses  a  fixed  Mcmd  for  the  desired  Mach  number.  AM  is  multiplied  by 
the  speed  of  sound  to  obtain  a  delta  velocity.  The  delta  velocity  is  divided  by  a  time 
constant  and  gravity  to  produce  a  dimensionless  gmach^  ■  Omachx  can  be  expressed  by 


Qmachx  —  i^Mcmd  M') 


(3.34) 


where 


a  =  the  speed  of  sound 
tm  =  the  time  constant 
g  =  gravity 


The  total  gcmdx  is  a  combination  of  gmachx  and  a  command  to  counter  the  effects  of 
gravity. 


Qy  ~  TeV 


0 

0 


-9\ 
—g  sin  7 

0 

L-5fcos7 


Making  gy^  dimensionless  and  opposite  to  the  gravity  term  yields 


gcmdgx  —  sm  7- 


(3.35) 


(3.36) 


(3.37) 
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The  total  command  in  the  horizontal  direction  is  expressed  by 

Qcmdx  —  Qmachx  d“  Qcmdgx 

=  (Mcmrf  -  M) -^  +  sin7.  (3.38) 

tmS 

ENGAGE  does  not  use  the  actual  value  of  the  speed  of  sound  in  air.  ENGAGE 
approximates  the  speed  of  sound  and  assumes  it  is  constant  with  a  =  304.3  m/sec 
(998.2  ft/sec). 

The  commands  in  the  vertical  direction,  z^,  are  a  combination  of  proportional 
navigation,  gravity  effect,  and  a  loft  command.  The  proportional  navigation  com¬ 
mand  is  expressed  in  Eq.  (3.30) 

V  Nj/ios 

9cmdP/N,  -  - - - • 

The  command  to  counteract  the  effect  of  gravity  in  the  Zy  direction  is  determined 
by  examining  the  gy  vector. 

■  0 

gy  =  Tev  0 

.-9 
—g  sin  7 

0  (3.40) 

—g  cos  7 . 

which  can  be  expressed  as  a  dimensionless  command  by  dividing  by  gravity 

i/cmdaz  =  cos  7.  (3.41) 

The  loft  command  is  employed  to  increase  the  altitude  of  the  missile,  giojt  is  the 
control  that  is  unknown  in  the  formulation  of  the  optimization  problem.  The  total 
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9cmdz  vector  is  expressed  as 


Qcmdz  —  9cmdP/Nz  d"  9cmdgz  d"  9loJt 

VN^ilos  , 

= - ]- COS 'f  +  gioff 

9 


(3.42) 


3.2  Control 

This  section  describes  the  control  logic  for  the  fuel  flow  rate  of  the  RAMJET 
engine,  the  missile  roll  control,  and  the  missile  pitch  control. 

3.2.1  Fuel  Flow  Rate  Control.  The  throttle  logic  accomplishes  three  things. 
First,  it  prevents  the  pressure  inside  the  combustion  chamber,  station  3  in  Fig.  2.11, 
from  getting  too  large.  In  turn,  this  prevents  too  large  a  pressure  at  the  inlet.  Ex¬ 
cessive  pressure  at  the  inlet  would  prevent  sufficient  air  flow  to  support  combustion. 
Second,  the  throttle  logic  maintains  Wf  between  a  minimum  and  maximum  value. 
Decreasing  Wf  too  low  will  cause  a  flame  out.  Increasing  Wf  too  much  will  result 
in  wasted  fuel  since  exceeding  the  stoichiometric  fuel  air  ratio  will  cause  fuel  to  be 
exhausted  without  being  converted  to  energy  by  the  combustion  process.  Third,  the 
throttle  logic  will  prevent  the  missile  from  exceeding  a  maximum  velocity.  For  sus¬ 
tainment  of  propulsion  in  the  combustor,  the  air  flow  rate  must  not  blow  the  flame 
out  the  rear  of  the  engine.  This  can  happen  if  the  velocity  of  the  air  through  the 
combustor  is  too  large. 

3.2.2  Roll  Control.  A  block  diagram  of  the  roll  control  for  the  missile  is 
pictured  in  Fig.  3.4.  The  control  logic  is  simply  a  feedback  of  the  actual  roll  angle, 
<7,  using  static  gains.  From  the  diagram, 

^  =  T-\^C^F„ 
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Figure  3.4  Block  Diagram  of  Roll  Control 


where 


a  =  the  roll  rate  of  the  missile, 
o'cmd  =  the  roll  command, 

Ca  =  the  feedforward  gain,  and, 
=  the  feedback  gain. 

Defining  as 


Eq.  (3.43)  is  reduced  to 

a  {t)  =  — llcTcmdllu-i  {t  -  to) - cr  {t)  (3.45) 

To  To 

where 


To 


Co 


1  +  FoCo 


(3.44) 


||<7'cmd||  =  the  magnitude  of  input 
u_i  =  the  unit  step  function 
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where 


u-i  {t  -  to) 


0,t  <  to 
l,t>  to 


Taking  the  Laplace  transform  of  both  sides  yields 


1  1 

sa  (s)  -  cr  (0)  +  —a  (s)  = - \\(Tcmd\ 

T„  S 


(3.46) 


Collecting  terms,  the  expression  becomes 


(s  +  — )  cr  (5)  =  - - llo-cmdll  +  CF  (0) 

\  Tff  /  T  (j  S 


(3.47) 


or,  in  terms  of  (t(s) 


,  1  e  ,  <^(0) 

<7  (sj  = - 7 - 7^  +  - 


(3.48) 


r,  +  3  +  i- 

The  first  term  on  the  right  side  of  the  equation  can  be  expanded  in  partial  fractions 
to  yield 

A  B  + 


■s  •s  +  “ 


T<t 


A  and  B  are  determined  from  the  above  expression  to  be 


(3.49) 


A  —  ||u’c7nd|| 

B  —  ll^cmdll' 


(3.50) 

(3.51) 


Substituting  these  values  for  A  and  B  yields 


||<7cmd||  __st„  ||<^cmd||  „-sto  , 


(3.52) 
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Figure  3.5  Roll  Response  to  a  5  Degree  Step  Command 
Inverting,  <t  as  a  function  of  time  becomes 

cr{t)  =  cr{0)e~^^  +  \\acmd\\u-i{t  -  Q  -  \\aomd\\e"^u-i{t  -  to).  (3.53) 

The  expression  for  cr{t)  is  continuous  for  all  time,  however  a{t)  is  discontinuous  at 
time  to-  Therefore,  the  commanded  value  of  roll  is  achieved  by  a  continuous  change 
in  a  rather  than  by  a  discrete  step.  The  first  two  terms  on  the  right  side  of  the 
expression  are  the  zero  state  solution  (7:7).  The  final  term  is  the  zero  input  solution 
(7:8)  and  assumes  the  response  is  due  to  initial  energy  in  the  system  and  the  input 
is  zero.  For  a  linear  system  as  the  roll  control  for  the  missile,  the  complete  response 
is  just  a  linear  superposition  of  the  zero  input  and  the  zero  state  response  (7:9).  The 
roll  response  for  a  step  input  of  5°  is  shown  in  Fig  3.5  assuming  the  input  occurs  at 
to  =  0  and  (To  =  0°. 

To  determine  (Jcmd,  the  command  vector  in  the  velocity  axis  is  transformed  to 
the  missile  body  axis  system. 

dcmds  '^VBdcmdy 
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Figure  3.6  Block  Diagram  of  Pitch  Control 


=  TvsTsBOcrndv 

'  9cmd^  COS  a  +  {gcmdy  sin  a  +  gcmd^  cos  (j)  sin  a 

—  9cmdy  COS  (J  Qcmdz  ^ 

.  -gcmda:  sin  Ot  +  {gcmdy  sin  <t  +  gcmd,  cos  cos  a  _ 
However,  for  a  BTT  missile  the  following  must  be  true: 


(3.54) 


gcmdy g  —  0  —  gcmdy  COS  (T  gcmdz  sin  O'. 


(3.55) 


The  expression  for  acmd  is  therefore 

(Tcmd  =  arctan  f  ]  .  (3.56) 

\gcmdz  J 

3.2.3  Pitch  Control.  The  pitch  control  of  the  VFDR  is  accomplished 
through  the  deflection  of  control  surfaces  on  the  tail.  The  deflections  cause  the 
rotation  of  the  body  to  a  new  angle-of-attack  and  lift  is  generated  by  the  wings 
and  the  body  (9:18).  The  angle  of  attack  of  the  missile  changes  proportionally  to 
the  difference  between  the  magnitude  of  the  commanded  acceleration  in  the  velocity 
yv  —  zv  plane  and  the  actual  acceleration  in  the  velocity  yv  —  zy  plane.  A  block 
diagram  of  the  pitch  control  is  pictured  in  Fig  3.6.  Therefore, 


d  —  ll^'cmdyzll  ^aeroyz 


(3.57) 
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The  magnitude  of  gcmdyz  is  given  by 


lIPcmdyzll  —  {dcmdy  "t"  9cmdz'} 


(3.58) 


and  the  magnitude  of  the  actual  acceleration  is  given  by 


n 


aeroyz | 


Fvy ( Fvz  ,  \ " 

w) 


•  n  V 

cos  a  sm  p  +  —  1  cos  a  > 

S 


ff-(^cosQsin/9+^)sintr} 


(3.59) 


Recalling  that  5  =  0  and  /3  =  0  for  the  VFDR,  the  equation  reduces  to 

T  .  _  L\^ 


\n 


aeroyz  j 


(sin"  ^  + cos V) 


4 


T  .  L 


(3.60) 


The  pitch  control  is  an  open  loop  system  and  integrates  the  command.  A  bounded 
input  does  not  produce  bounded  output.  An  angle  of  attack  response  for  a  Ig  delta 
command  is  pictured  in  Fig  3.7.  As  is  the  case  with  the  roll  control  logic,  the  pitch 
control  logic  produces  a  continuous  change  in  the  angle  of  attack  so  discontinuities 
in  a  do  not  exist. 
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Figure  3.7  Angle  of  Attack  Response  to  Ig  Command 
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IV.  Analysis  of  Simulation 

Upon  completion  of  the  modeling,  the  simulation  was  coded  into  MATLAB^^. 
Two  baseline  test  flights  were  used  for  comparison  and  the  results  of  ENGAGE  were 
not  duplicated  for  flight  based  solely  on  proportional  navigation.  For  the  first  test 
flight,  the  missile  and  target  are  initially  at  the  same  altitude  and  relative  cross 
range;  the  flight  is  in  a  vertical  plane.  Figure  4.1  depicts  the  scenario.  The  target 
flies  a  constant  altitude  profile  and  is  non  accelerating  at  an  air  speed  of  0.9  Mach. 
The  missile  is  launched  with  an  initial  speed  of  0.9  Mach  and  is  pointing  directly  at 
the  target.  For  the  second  test  flight,  the  target  is  at  a  lower  altitude  and  a  different 
cross  range  position  than  the  missile.  Table  4.1  summarizes  the  initial  conditions  and 
the  results  for  the  two  flights.  The  states  for  the  first  flight  are  shown  in  Figs  4.2 
and  4.3.  The  solid  lines  indicate  the  simulation  results  from  the  MATLAB™  code 
and  the  dashed  lines  represent  the  results  obtained  from  ENGAGE.  An  error  exists 
between  the  two  simulations.  The  states  for  the  second  test  case  are  shown  in  Figs  4.4 
and  4.5.  Again,  an  error  is  noticeable. 

During  the  course  of  comparing  results  to  those  obtained  from  ENGAGE,  two 
errors  were  found  in  the  original  code  of  ENGAGE.  The  first  error  was  in  the  model 
of  the  atmosphere.  Recall  the  two  expressions  obtained  for  density  depending  upon 


Figure  4.1  Initial  Conditions  of  Baseline  Simulation 


Flight  1 


Flight  2 


Missile 

Target 

Missile 

Target 

Down  Range  (km) 

0 

27.8 

0 

18.5 

Cross  Range  (km) 

0 

0 

0 

9.26 

Altitude  (m) 

6098 

6098 

6098 

4573 

Mach  Number 

0.9 

0.9 

0.9 

0.9 

Flight  Path  Angle  (Deg) 

0 

0 

0 

0 

Heading  Angle  (Deg) 

0 

-180 

0 

-180 

Intercept  Time  (Sec) 

24.76 

21.10 

Rlos  (ni) 

2.03 

0.640 

Table  4.1  Initial  Conditions  and  Results  for  Baseline  Simulations 


Angle  of  Attack  (Deg) 


Time  (sec)  Time  (sec) 


Time  (sec)  Time  (sec) 

Figure  4.3  7,  ip,  a  and  cr  Simulation  Comparisons  for  Flight  1 
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Angle  of 


—  Original  Density  Profile : 
-i-  -  Corrected  Density  Profile 


0.2  0.4  0.6  0.8  1  1.2 

Density,  [kg/cu.  m] 

Figure  4.6  Density  as  a  Function  of  Altitude 
the  characteristics  of  temperature  for  the  local  region;  Eq.  (2.15)  and  Eq.  (2.16). 

/r„y(‘+*)  dT 

po  =  Pb\JfJ  for  ^  T  and, 

8(>»g-H)  dT 

Po  =  Pbf^  for  =  0. 

dha 


The  first  expression  is  used  to  approximate  air  density  in  regions  that  temperature 
is  characterized  as  varying  linearly  with  altitude,  Fig  2.1.  The  second  expression  is 
used  within  regions  that  temperature  is  constant  with  altitude. 


ENGAGE  utilizes  the  expression  for  density  valid  only  in  isothermal  regions. 
Figure  4.6  shows  the  density  calculated  with  both  equations.  Figure  4.7  shows  the 
percent  error  between  the  correct  expression  and  the  incorrect  expression.  Notice 
that  near  the  top  of  the  gradient  region  the  error  is  almost  10.0  percent.  The  density 
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2  3  4  5  6  7  8 

Percent  Error 


Figure  4.7  Percent  Error  in  the  Density  Value 


value  is  used  to  determine  the  aerodynamic  forces  acting  on  the  missile;  Eqs.  (2.54) 
and  (2.57)  where 

<3.  =  \pV\  (4.1) 

Hence,  the  aerodynamic  forces  acting  on  the  missile  in  a  region  of  the  atmosphere 
that  temperature  is  modeled  as  varying  with  altitude  will  be  in  error  and  this  error 
could  be  a  high  as  9.7  percent  at  8,400  m. 

The  throttle  control  logic  during  the  RAMJET  boost  phase  requires  the  den¬ 
sity  of  the  air  in  determining  the  maximum  airflow  rate  through  the  engine.  The 
maximum  airflow  rate  information  is  then  utilized  in  the  determination  of  the  fu- 
elflow  rate.  Consequently,  an  incorrect  density  value  will  result  in  an  improper  thrust 
level. 

The  thrust  calculated  during  the  rocket  boost  phase  will  also  have  an  associated 
error  due  to  the  density  error.  Recall  Eq.  (2.59),  the  thrust  for  the  rocket  boost  phase. 


Ft  —  Tgl  -f-  Aexiti^Psl  p}‘ 


(4.2) 


Due  to  the  error  in  the  density  calculation,  the  atmospheric  pressure  calculated  will 
be  lower  than  the  actual  atmospheric  pressure.  Hence,  the  thrust  value  calculated 
in  ENGAGE  will  be  too  large. 

A  second  error  was  discovered  in  the  line  of  sight  algorithm  of  ENGAGE. 
Recall  from  Eq.  (3.21)  that 


7LOS  = 


1  {^los  (®LOS  +  ylos) 


Flos i^los  +  vios)^ 

-ZLOS  {yxLos^LOS  +  ^VLOsyLOs)}  ■ 
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Figure  4.8  Percent  Error  in  'yios 
ENGAGE  calculated  ^os  as 


7los  = 


^LOS  i^LOS  +  VlOS  +  ‘^^LOsY 
-ZlOS  {V^los^LOS  +  VyLosyLOs)]  ■ 


T  {^ZLos  {^LOS  +  VlOs) 


Upon  inspection,  an  error  between  the  two  expressions  for  jias  occurs  when  the 
target  and  missile  are  not  at  the  same  altitude.  The  relative  error  between  the  two 
expressions  is  shown  in  Fig  4.8.  The  x^os  and  y^os  values  were  both  fixed  at  9.3  km 
and  zios  was  varied  from  -8  to  8  km.  Recall  the  proportional  navigation  command 
in  the  vertical  direction: 


V  NjfLos 

dcmdP/N^  =  - - - 
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The  command  is  proportional  to  jLOS-  Hence,  ENGAGE  generated  less  of  a  vertical 
proportional  navigation  command.  Figure  4.8  depicts  the  percentage  decrease  of  the 
proportional  command. 

The  errors  noticed  in  Figs  4.2  through  4.5  are  a  combination  of  these  two 
errors.  Modifying  the  MATLAB^^  code  to  recreate  the  errors,  the  simulations 
agreed  and  the  missile  still  intercepted  the  target.  The  MATLAB™  simulation  can 
interface  with  computer  workstations  so  that  numerical  techniques  can  be  employed 
to  optimize  the  missile’s  performance. 
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V.  The  Optimal  Control  Problem 

5.1  Necessary  Conditions  For  a  Minimum 

This  chapter  describes  the  necessary  conditions  for  the  standard  optimal  con¬ 
trol  problem  with  free  final  time  (2:71)  (5:155).  Consider  a  dynamic  system  described 
by  a  set  of  differential  equations  with  known  initial  time,  to,  in  the  form  of 

X  =  f{x,u,t)  (5.1) 

where, 

X  =  an  n  X  1  vector  of  the  state  variables, 
u  =  an  m  X  1  vector  of  control  variables,  and 
to,Xo=  known. 

The  system  is  subject  to  p  constraints  on  the  terminal  conditions  of  the  state  vari¬ 
ables,  q  constraints  on  the  state  variables,  and  r  constraints  on  the  control  variables, 
so  the  constraints  can  be  expressed  by 


^,(a;/,t/)  =  0  for  i  =  1  . .  .p,  (5.2) 

5j(a:,f)  =  0  for  f  =  1  . .  .9,  and  (5.3) 

Ci{u,  f )  =  0  for  f  =  1  . . .  r.  (5.4) 

The  problem  is  to  determine  a  function  u{t)  which  minimizes  a  scalar  performance 
index,  J,  subject  to  the  conditions  stated  in  Eqs  (5.1)  through  Eq.  (5.4).  The 
performance  index  is  given  by 

J  =  ^{xf,tf) -\-  f  ^  L{x,u,t)dt.  (5.5) 
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The  performance  index  is  adjoined  with  the  system  differential  equations,  Eq.  (5.1), 
by  introducing  Lagrange  multipliers,  A(t),  as  a  multiplier  function.  Similarly,  the 
constraints,  Eqs  (5.2)  through  (5.4),  are  also  adjoined  to  the  performance  index  with 
additional  Lagrange  multipliers,  d(t),  and  Therefore,  the  augmented 

performance  index,  J',  is  given  by, 


J'  =  G{tf,Xf,u)  + 


X 


dt 


(5.6) 


where  the  Bolza  function,  G,  and  the  variational  Hamiltonian,  H  (5:156),  are  defined 
as 


Xf,u)  =  $(t/,  Xf)  +  i/^^(t/,  Xf), 

H{t,x,u,X,0,fi)  =  L{t,x,u)  +  X'^f{x,t,u)  +  0’^S{x,t)  +  fjFC{u,t). 

Furthermore,  the  constraints  can  be  of  an  equality  or  inequality  type.  For  inequality 
constraints,  slack  variables,  /,  are  introduced.  For  example,  assume  a  control  u  is 
limited  to  some  maximum  value,  Umax-  The  constraint  equation  is, 

U  ^  Umax-  (^*^) 

The  new  constraint  equation  is  expressed  by  Eq.  (5.8), 


U  -  Umax  +  1^  =  0  (5.8) 

The  inequality  constraint  is  now  an  equality  constraint  because  any  value  of  u  greater 
than  Umax  will  never  satisfy  Eq.  (5.8)  for  real  valued  1.  When  I  is  zero,  the  control 
is  at  the  constraint  value.  When  I  is  not  zero,  is  identically  zero.  A  potential 
drawback  is  there  exists  another  unknown  variable.  Note  the  Hamiltonian  is  now  a 
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function  of  I  also.  Thus, 

H{t,x,  u,X,6,  /)  =  L{t,  X,  u)  +  f{x,  t,  u)  +  6^S{x,t)  +  ^^C{u,  t,  /)(5.9) 

and  the  augmented  performance  index  becomes 


J' =  G{tf,Xf,if)  +  {t,x,u,X,0,  fx,l)  —  X^x  dt.  (5.10) 


The  first  variation  of  the  augmented  performance  index  is  given  by 


6J'  —  GijStj  +  GxjSxf  +  Gi,8v 

+  C\Hx8x  +  hJu  +  Hx  ~SX  +  He60  +  hJii  +  Hi  ~8l  -  -  X^8x]dt 

Jto 

+[H  -  A’-ilft  \\i  (5.11) 


where  the  notation  implies 


()v  = 


dj) 

dy  ’ 


o/\ 

{)tj  =  evaluated  ai  t  =  tj, 


8{)  =  a  time  free  variation, 
^0  =  a  time  fixed  variation. 


The  last  term  in  Eq.  (5.11)  can  be  evaluated  at  to  and  tf  as 

[H  -  X^x]8t  \li=  {Hf  ~  X'^jXf)8tj  -  {Ho  -  Xlxo)8to.  (5.12) 

The  initial  time  is  known,  hence  8to  is  zero.  Therefore,  the  final  term  on  the  right 
in  Eq.  (5.12)  is  eliminated.  Eliminating  this  term,  the  variation  of  the  augmented 
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performance  index  becomes 


SJ'  =  GtjStf  -f  GxjSxf  +  GySv  H-  {Hf  -  X'^Xf)Stf 

+  f'^HxSx  +  hJu  +  Hx8X  +  HeSe  +  H^Sti  +  Hi~8l  - 
Jio 

x^8\  —  X'^6x]dt. 

(5.13) 

Recall  the  constraint  equations  given  by  Eqs  (5.2)  through  (5.4). 

Recognizing 

G,  =  ®(®/,fy)=0 

(5.14) 

He  =  S{x,t)  =  0 

(5,16) 

i/^  =  C(«,t)  =  0 

(5.16) 

eliminates  HgSO  and  from  Eq.  (5.13).  Thus,  the  variation  of  the 

augmented  performance  index  is  now  expressed  as 

SJ'  =  GtfStf  +  GxfSxf  +  (Hj  —  X^Xf)6tf 

+  f'  [Hx'Sx  +  Hju  +  HiU  +  {Hx  -  x^)SX  -  X^lx]  dt.  (5.17) 

V  tc 

Noting  Hx  =  and  recalling  Eq.  (5.1),  the  terms  multiplied  by  reduce  to 

Hx  —  x^  =  —  x^  =  x^  —  x^  =  0.  (5.18) 

Hence,  the  term  multiplied  by  6X  is  always  zero.  Substituting  Eq.  (5.18)  into 
Eq.  (5.17)  yields 

6f  =  Gt^6tf  +  Gxf6xf  +  {Hj-X'^jif)6tf 

+  {HxSx  +  hJu  +  Hi~6l  -  X'^6x}dt.  (5.19) 
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X 


Nominal  Path 
Neighboring  Path 


to  tftf  +dt^ 

Figure  5.1  Relationship  Between  6xf,  8xf,  and  XfSt 

The  last  term  in  the  integral  in  Eq.  (5.11)  can  be  integrated  by  parts  in  the  following 


fashion, 


/  ^  —\^6x  =  —X'^Sxf  +  XgSxo  +  f  ^  —X^Sxdt. 
Jto  ■'  Jto 


(5.20) 


Therefore,  substituting  the  above  information  into  Eq.  (5.11),  6J'  can  now  be  ex¬ 
pressed  by 

SJ'  =  GtjStf  +  GxjSxf  -  XjSxf  -h  X^Sxo 

+  {Hf  -  XjxAStf  +  [(H^  +  X^)Sx  +  HJu  +  H,Sl]dt.  (5.21) 

Jtn 


The  relationship  between  the  final  time  free  variation  and  the  final  time  fixed  varia¬ 


tion  is 


8xj  =  8xf  —  XfSt. 


(5.22) 


This  relationship  is  presented  graphically  in  Fig  5.1.  Also,  since  the  initial  conditions 


are  known. 


Sxo  =  0  X^SXo  =  0. 


(5.23) 
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Substituting  Eq.  (5.22)  and  Eq.  (5.23)  into  Eq.  (5.21)  yields 


SJ'  -  Gt^Stj  +  Ga;^6xf  +  X^XfStj  -  \'^6xf  +  [Hj  -  \'^Xf]Stf 
+  f  ^  [{H^  +  X^)Sx  +  Hju  +  Hi6l]dt. 

to 

Collecting  terms  and  noting  Aji  -  AJ i  =  0  yields 

6J'  =  [Gtj  +  Hf)Stf  +  (Gx^  —  X'j)6xf 

+  j  '  [{Hx  +  X^)Sx  +  HJu  +  Hi8l]dt.  (5.24) 

Jto 

For  an  extremal,  the  variation  of  J'  must  be  zero.  Additionally,  each  variation  is 
independent  yielding  the  free  final  time  necessary  conditions  for  an  extremal; 


1 

II 

(5.25) 

> 

II 

(5.26) 

^  =  -Hl 

(5.27) 

iT«  =  0 

(5.28) 

JT,  =  0 

(5.29) 

^  =  S  =  C  =  0 

(5.30) 

^0  —  ^05 

(5.31) 

—  ^os 

(5.32) 

X  =  f, 

(5.33) 

Equations  (5.25)  through  (5.33)  provide  the  necessary  conditions  for  an  extremal. 

The  Weierstrass  condition  and  the  Legendre-Clebsch  condition  are  two  suffi¬ 
cient  conditions  for  a  minimum  (5:106).  The  Weierstrass  condition  is  necessary  for 
a  strong  relative  minimum.  This  condition  states  that  no  other  control  can  pro¬ 
duce  a  lower  value  of  the  Hamiltonian  than  the  optimal  control.  More  precisely,  the 
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Weierstrass  condition  states  that  for  any  control  u*  other  than  the  optimal  control 
u 

H(t,x,U:t,X)  —  H{t,x,u,X)  >  0.  (5.34) 

The  Legendre- Clebsch  condition  is  sufficient  for  only  a  weak  relative  minimum. 
The  weak  variation  is  a  particular  case  of  a  strong  variation  and  is  applicable  only 
for  small  variations  in  the  control.  The  Legendre-Clebsch  condition  requires 

H^u  >  0.  (5.35) 

When  the  Weierstrass  condition  is  satisfied,  the  Legendre-Clebsch  condition  is  auto¬ 
matically  satisfied.  However,  the  Weierstrass  condition  is  not  necessarily  satisfied  if 
the  Legendre-Clebsch  condition  is  satisfied.  When  the  Legendre-Clebsch  condition 
is  not  satisfied,  the  Weierstrass  condition  will  not  be  satisfied. 

5.2  Shooting  Method 

The  shooting  method  is  an  algorithm  that  provides  a  numerical  solution  to  a 
two-point  boundary-value  problem.  Hence,  if  the  optimal  control  problem  can  be 
converted  into  a  two-point  boundary-value  problem,  the  shooting  method  can  be 
used  to  provide  a  numerical  solution.  Recall  the  necessary  conditions  for  the  free 
final  time  problem; 


x  =  f 

(5.36) 

(5.37) 

^0  -  ^OS 

(5.38) 

Xo  -  Xos 

(5.39) 

Hf  = 

(5.40) 
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0  =  Hu 

(5.41) 

A/  = 

(5.42) 

Hi  =  0 

(5.43) 

^  =  S  =  C  =  0. 

(5.44) 

The  Lagrange  multipliers  adjoined  to  the  constraints,  i/,  B  and  ^  and  the  slack 
variables,  i,  are  eliminated  through  the  solutions  to  Eqs  (5.40)  through  (5.44).  The 
control,  w,  is  expressed  as  a  function  of  «,  A,  and  t  from  the  solution  to  Eq.  (5.41). 
Hence  the  problem  is  reduced  to  a  function  of  a;.  A,  and  t.  The  necessary  conditions 
become 


x  =  Ft  {t,x,X) 

(5.45) 

X  =  F2  {t,x,X) 

(5.46) 

h{tf,Xf,Xf)  =  0 

(5.47) 

io  ^os 

(5.48) 

Xo  =  Xos 

(5.49) 

A  new  vector  is  defined  that  is  a  combination  of  x  and  A 

r«i 


z  = 


(5.50) 


Therefore,  the  necessary  conditions  are  restated  in  terms  of  this  new  definition  as 


z  =  F{t,z) 


io  —  ^os 


(5.51) 

(5.52) 


®o 


— 


Ao 


h{tf,Zf)  =  0. 


(5.53) 

(5.54) 
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z  is  a  vector  of  length  2n  where  n  is  the  number  of  original  states,  n  of  the  initial 
conditions  of  z  are  known  (i.e.  Xo).  Equation  (5.54)  gives  the  2n  +  1  terminal 
conditions  for  Zf  and  tf.  The  conditions  stated  in  Eqs  (5.51)  through  (5.54)  represent 
a  standard  two-point  boundary-value  problem.  The  n  unknown  initial  z  values  and 
a  final  time  are  guessed  and  the  shooting  method  iteratively  updates  the  unknown 
initial  z  (i.e.  Ao)  and  until  Eq.  (5.54)  is  satisfied  to  within  a  specified  tolerance. 
Therefore,  a  method  is  required  to  update  the  initial  guesses. 

Taking  the  first  variation  of  Eq.  (5.51)  yields 

8z  =  FJz  (5.55) 

J 

-j-Sz  —  F^Sz.  (5.56) 


Recall  the  definition  for  z  and  that  the  initial  conditions  for  x  are  known.  Therefore, 
the  total  variation  of  z  at  the  initial  time  is 


Sxo 

SXo 

0 

8\o 


(5.57) 

(5.58) 


Assume  the  initial  guesses  are  in  the  neighborhood  of  the  actual  values.  A  state 
transition  matrix,  $,  is  introduced  to  relate  the  total  variation  of  z  at  any  time  to 
the  initial  total  variation  of  z.  Using  $,  6z  is  expressed  by 


^z  =  $  {t,to)SZo 


(5.59) 


and  the  initial  condition  on  $  is 


^  =  I 


(5.60) 
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Substituting  Eq.  (5.59)  into  Eq.  (5.56)  yields 


—  [$  {t,  to)  Szo\  ^  (i,  to)  Szo 

$  (t,  to)  Szo  =  (t,  to)  Szo. 


(5.61) 

(5.62) 


Equation  (5.62)  reduces  to  another  set  of  differential  equations 


$  (tjo)  =  F^^  (t,to) 


(5.63) 


with  the  initial  condition 


^{to,to)  =  J. 


(5.64) 


This  equation  is  integrated  from  to  to  tj  to  yield  an  expression  for  the  final  total 
variation  of  z  in  terms  of  linear  combinations  of  the  initial  total  variation  of  z 


Szj  —  $  (^/i^o)  ^^0 


(5.65) 


Applying  the  initial  total  variation  of  z  to  the  above  equation,  the  final  total  variation 
of  z  becomes 


6zf  =  ^  {tf,to)6zo 


=  ^{tf,io) 


(5.66) 


The  transition  matrix  can  be  partitioned  to  yield 


hf  =  [^l{tf,to)  ^2{tf,to)] 


(5.67) 
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Therefore,  the  final  total  variation  of  z  is  expressed  solely  as  a  total  variation  of  the 
initial  unknown  elements  of  z 

62/ =  $2  (5.68) 

The  final  boundary  conditions  of  Eq.  (5.54)  are  evaluated  based  upon  the 
initial  guess  for  the  unknown  initial  variables.  Assuming  Eq.  (5.54)  is  not  satisfied, 
the  initial  guess  of  the  unknown  variables  must  be  modified.  To  determine  the 
modification  to  the  guess,  the  time  free  variation  of  the  terminal  boundary  condition 
is  examined 

6h  =  +  hzjSzf.  (5.69) 

The  relationship  between  a  time  free  variation  and  a  total  variation  is 

6zf  =  6zf  +  ZfStf.  (5.70) 

Substituting  Eq.  (5.70)  into  Eq.  (5.69)  yields 

Sh  =  htjStf  +  h^^Szf  +  h^^ZfStf 

=  h^jSzf  +  (htj  +  h^jZj^  Stj.  (5.71) 

The  total  variation  of  the  final  z  is  known  through  Eq.  (5.68).  Making  the  substi¬ 
tution  for  Szf,  Sh  becomes 

Sh=  h^^^2{tf,io)SXo+  (ht^  +  Stf.  (5.72) 

The  variation  of  h  can  also  be  expressed  as 

Sh  =  hnew  —  h.  (5.73) 
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Recall  the  previous  assumption  of  the  neighboring  path.  With  this  assumption,  the 
variation  of  h  due  to  variations  in  z  is  small.  Hence,  the  value  for  hne^i  is  close  to 
the  old  value  of  h  and  is  expressed  as 

hnew  =  rjh  (5-74) 

with  T)  close  to  1.  (5.75) 

Substituting  the  expression  for  knew  into  Eq.  (5.73),  6h  becomes 


8h  =  —  {I  —  rj)  h 


=  —ah 


(5.76) 


where  a  is  defined  as. 


a  =  (1  -7?). 


(5.77) 


Typically,  a  is  about  0.1.  Combining  Eq.  (5.72)  and  Eq.  (5.76)  yields 


[  kzj^2  (f/i  to) 


ktj  4"  hzj Zj) 


SXo 

6tf 


—ah. 


(5,78) 


A  means  to  modify  the  initial  guess  is  provided  by  Eq.  (5.78).  The  initial  guess  for 
Ao  and  f/  is  updated  by 


Anew  —  Ao  +  S\o  (5.79) 

=t}  +  Stf  (5.80) 

To  implement  the  algorithm,  Ao  and  tf  are  guessed.  Equations  (5.51)  and  (5.63)  are 
integrated  from  to  to  tf.  Zf,  $2(^/5 to),  and  \\h\\  are  computed.  \\h\\  is  checked  to 
determine  if  it  is  within  a  specified  tolerance,  a  is  set  to  1  and  ^Ao  and  St j  are  cal¬ 
culated  and  intermediate  values  for  the  initial  guesses  are  formed.  Equations  (5.51) 
and  (5.63)  are  integrated  again  from  to  to  tf  using  the  intermediate  values.  \\h\\  is 
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Figure  5.2  Lunar  Launch  Example  Problem 


checked  again  to  determine  if  its  value  decreased.  If  not,  a  is  reduced  by  half  and 
SAo  and  Stf  are  recomputed.  If  \\h\\  was  less  than  the  previous  value,  then  the  inter¬ 
mediate  values  are  accepted  as  the  new  Ao  and  tj  values.  This  process  is  repeated 
until  the  specified  tolerance  for  \\h\\  is  achieved. 

5.3  Shooting  Method  Example  Problem 

The  necessary  conditions  are  utilized  to  formulate  a  two-point  boundary- value 
problem  for  a  lunar  launch.  The  shooting  method  previously  described  is  imple¬ 
mented  to  provide  a  numerical  solution  to  the  problem.  The  problem  is  to  find  the 
control  history  which  minimizes  the  final  time  for  a  spacecraft’s  injection  into  Lunar 
orbit  at  prescribed  final  conditions  from  the  surface  of  the  moon.  Figure  5.2  depicts 
the  scenario.  The  problem  is  stated  as; 


minimize  J  —  tj 


(5.81) 


5-13 


subject  to  the  differential  equations 


X  —  u 

(5.82) 

y  =  v 

(5.83) 

T  ^ 
u  =  —  cos  6 

(5.84) 

m 

m 

T 

V  =  —  sin  9  —  g 
m 

(5.85) 

with  initial  conditions 


^0  =  0 

(5.86) 

a:(0)  =  0 

(5.87) 

y(0)  =  0 

(5.88) 

«(0)  =  0 

(5.89) 

v{Q)  =  0 

(5.90) 

and  the  terminal  conditions 


tf  =  unknown 

(5.91) 

x{tf)  =  free 

(5.92) 

y{tf)  =  15.24  km 

(5.93) 

u{tj)  =  1660  m/s 

(5.94) 

u  {tj)  =  0 

(5.95) 

where 


u  =  the  velocity  component  in  the  x  direction 
V  =  the  velocity  component  in  the  y  direction 
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T 

—  =  6.34  m/s^  and  is  constant 
m 

g  =  the  force  of  gravity,  1.62  m/s^  and  is  constant 
6  =  the  control  variable. 


Expressing  the  terminal  constraints  as  functions  of  ^  yields 


=  15,240  =  0 

(5.96) 

=  u/  —  1660  =  0 

(5.97) 

c» 

II 

II 

p 

(5.98) 

Adjoining  the  terminal  constraints  with  Lagrange  multipliers,  i/,  the  Bolza  function 
is  expressed  as 

G  {tf,  Xf,  u)  =  <t)  {tf,  Xf)  +  {tf,  Xf) 

=  tf  +  t^y{yf-  15240)  +  -  1660)  +  i/vVf.  (5.99) 

Adjoining  the  differential  equations  with  Lagrange  multipliers.  A,  the  Hamiltonian 
is  expressed  as 


H  (t,  X,  0,X)  =  L  {t,x,0)  +  A^/  (<,  X,  6) 

=  +  A„  ( —  cos  +  A„  ( —  sin  0  — 

\m  /  \m  / 


(5.100) 


The  augmented  performance  index  is  given  by 


J' =  G{tf,Xf,v)  +  {t,x,u,\)  —  X^x  dt 

=  tj-{-Uy{yf  -  15240)  +  i/„  (u/  -  1660)  +  VyVj 

+  ’  j^Aj;U  +  \yV  +  A„  cos  ^ 

+  Ai,  f  —  sin  ^  —  o')  —  A^x  dt. 

\m  J 


(5.101) 
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The  necessary  conditions  for  a  minimum  result  from  the  first  variation  of  the  aug¬ 
mented  performance  index  being  set  equal  to  zero.  Hence,  the  differential  equations 
for  A  result  from  the  necessary  condition 


(5.102) 

Applying  Eq  (5.102),  the  expressions  for  A  are 

o 

II 

•-< 

(5.103) 

Ay  —  0 

(5.104) 

<< 

1 

II 

(5.105) 

A^)  —  Aj^. 

(5.106) 

The  control  is  eliminated  by  the  necessary  condition 

II 

p 

(5.107) 

Applying  Eq.  (5.107)  yields 

Hg  =  —Xu  ( —  sin^')  +  Ay  ( —  cos 

\m  )  \m  J 

(5.108) 

=  0. 

(5.109) 

Solving  for  9  yields  two  possible  control  laws 

A„  sin  9  =  Xy  cos  9 

(5.110) 

sm9  Xy 

=  ia,n9  = 
cos  9  Xu 

(5.111) 

/I 

ui  =  arctan  — 

Art 

or  92  =  arctan  — — 

A^^ 

(5.112) 

(5.113) 
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I 


I 


The  quadrant  problem  shown  in  Fig  5.3  results  because  of  the  two  possible  solutions 
for  9.  Evaluating  the  second  partial  derivative  of  the  Hamiltonian  with  respect  to  6 
resolves  which  solution  is  sufficient  for  a  local  minimum.  Taking  the  second  partial 
derivative  of  H  with  respect  to  6  yields 


d^H 

W 


T  T 

=  —  A„ —  cos  9  —  Xy —  sin0 
m  m 

-T 

= - (A„  cos  9  Xy  sin  9) . 

m 


Substituting  9\  into  Eq.  (5.114)  yields 

-T  (  XI  ^  a;  \ 

aip  m  \(a;  +  a;)5  (aj  +  a;)*) 

<  0. 


(5.114) 


(5.115) 

(5.116) 
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di  is  a  local  maximum.  Substituting  $2  into  Eq.  (5.114)  yields 

_-T  !  -XI  XI  \ 

“  l(A5  +  AJ)i  (K  +  yj) 

>  0. 


(5.117) 

(5.118) 


$2  satisfies  the  sufficient  condition  for  a  local  minimum  because  it  is  always  positive. 
Hence,  6  is  eliminated  from  all  equations  by  using  the  relationships 


sin0  = 


005  02  = - — — r- 

(K  + 


The  terminal  conditions  are  the  three  constraint  equations  and  =  0  because 
x{tf)  is  unspecified.  The  necessary  condition  Hj  =  -Gt^  provides  the  additional 
expression  required  to  determine  t /.  These  conditions  are  expressed  by 


{"t  f)  U  {t  f)  \y  {t  f)  V  {t  f)  g\y[tj)  +  l 

O' 

y(t/)/15, 240-1 

0 

u(f/)/1660  -l 

= 

0 

V  (tf)  /1660 

0 

K  (tj) 

.0. 

where 

A  =  (a2  +  A2)"  (5.119) 

The  problem  is  now  in  the  form  of  a  standard  two-point  boundary- value  prob¬ 
lem.  To  obtain  a  numerical  solution  to  the  problem  using  the  shooting  method,  the 
lunar  launch  problem  is  formulated  to  coincide  with  the  description  of  the  shooting 
method  provided  in  the  previous  section.  Defining 


z  =  [x  y  u  V  \y 


(5.120) 
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the  differential  equations  become 


z  = 


u 


m  A 

-  — ^  —  0 
m  A  “ 

0 

0 

'^x 

A,, 


(5.121) 


where  A  is  defined  by  Eq.  (5.119).  The  initial  conditions  are  given  by 


=  0 


rp 

z((,)  =  10  0  0  0  A  xo  ^yo  A^o  A^^  j  . 


(5.122) 

(5.123) 


The  terminal  conditions  are  expressed  in  the  form  of  Eq.  (5.54)  by 


■  X.  (t,) «((,)  +  A,  +  r 

»((/)/15240  -l 
u  (tf)  /1660  —  1 
u(t/)/1660 
A,(</) 


Note  some  scaling  has  been  applied  to  the  terminal  conditions  for  numerical  condi¬ 
tioning.  The  differential  equations  for  the  state  transition  matrix,  Eq.  (5.63),  require 
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an  expression  for  Fz  which  is  expressed  as 


where 


Fz  = 


0  0  10 
0  0  0  1 
0  0  0  0 
0  0  0  0 
0  0  0  0 
0  0  0  0 
0  0  0  0 
0  0  0  0 


T  (-1  XV 

m  I  A  A3 , 


B  = 


T  Xu^v 


m  A3 

m 


A  A3. 


0 

0 

0 

0 

0 

0 

-1 

0 


0  0  O' 

0  0  0 
0  A  B 
0  B  C 
0  0  0 
0  0  0 
0  0  0 
-10  0 


(5.124) 


(5.125) 

(5.126) 

(5.127) 


The  update  to  the  initial  guesses  requires  expressions  for  ht^  and  hzj  and  these  are 
expressed  by 


htf  =  [0  0  0  0  of 


hzj  — 


where  D  = 

E  = 


(5.128) 


'0 

0 

Xx{tj) 

Ay(ty)  U 

(tf) 

V{tf) 

D{tj) 

EitfY 

0 

1 

15240 

0 

0 

0 

0 

0 

0 

0 

0 

1 

1660 

0 

0 

0 

0 

0 

(5.129) 

0 

0 

0 

1 

1660 

0 

0 

0 

0 

.0 

0 

0 

0 

1 

0 

0 

0  . 

T 

m 

A=> 

t>l 

1 

2A„ 

A 

(5.130) 

T 

m 

'KK  ^ 

A’  + 

A3 

2A„ 

A 

-9- 

(5.131) 

Applying  the  shooting  method  to  the  sample  lunar  launch  problem,  the  min¬ 
imum  final  time  was  determined  to  be  272.706  seconds.  The  tolerance  for  ||/i||  was 
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Iteration 


A 


Table  5.1  Iteration  History  of  Ao  for  Shooting  Method  Example  Problem 


Iteration 

tf  (Sec) 

Initial 

300 

1 

240.08 

2 

271.45 

RfciirariBi 

3 

272.68 

4 

272.71 

5 

272.71 

Table  5.2  Iteration  History  of  and  \\h\\  for  Shooting  Method  Example  Problem 

1  X  10“®  and  was  achieved  in  five  iterations.  Table  5.1  summarizes  the  iteration 
history  of  Aq.  Table  5.2  summarizes  the  iteration  history  of  tf  and  \\h\\.  The  time 
histories  of  the  states  are  shown  in  Fig  5.4.  Table  5.3  summarizes  the  prescribed 
final  states  and  the  numerical  results  for  the  final  states.  The  time  history  of  the 
control  required  to  achieve  the  terminal  conditions  is  shown  in  Fig  5.5. 


®(^/) 

Prescribed  Value 

Numerical  Result 

x{tf) 

Free 

222.29723  (km) 

y{tj) 

15.24  (km) 

15.239999  (km) 

U{tf) 

1660  (m/s) 

1660.0000  (m/s) 

V{tf) 

0 

-8.86199x10-^  (m/s) 

Table  5.3  Summary  of  x{tf)  for  Shooting  Method  Example  Problem 
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Time  (Sec) 


Time  (Sec) 


Figure  5.5  Control  History  for  Lunar  Launch 


5.4  Three-Dimensional  Pursuit 


Consider  the  three-dimensional  intercept  scenario  where  the  problem  is  to  find 
the  control  u{t)  that  minimizes  the  final  time.  The  performance  index  is  given  by 


J  —  tf. 

The  differential  equations  described  in  chapters  2  and  3  are 


(5.132) 


where 


and  the  terminal  constraint 


X  =  V  cos  7  cos  V’ 

(5.133) 

y  =  V  cos  7  sin  ^ 

(5.134) 

h  =  V  sin7 

(5.135) 

V  =  cos  a  —  D)  —  g  sin  7 

(5.136) 

7  =  (T  sin  q;  H-  T)  cos  <T  ^  cos  7 

(5.137) 

Ip  =  - (T  sin  a  +  I)  sin  a 

VW  cos  7 

(5.138) 

a  =  C^  -f.  ^  (T  sm  a  +  L) 

(5.139) 

If,  f  9cmdy 

<T  =  —  arctan  ( - 1  —  <t 

L  \  u  / 

(5.140) 

W  =  —Wf 

(5.141) 

[  l,a>0 
phase =  < 

(  -l,a  <  0 

(5.142) 

V>(t/)=  {xf~XTjy+[yf-yT,Y-\-(hf-hTjy'^=0.  (5.143) 
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The  control,  n,  is  chosen  to  represent  the  total  commanded  g’s  in  the  zy  direction. 
Recalling  Eq.  (3.42),  u  is  defined  as 

^  —  Qcmdz  —  9cmdP/Nz  T  Qcmdgz  T  9loft‘  (5.144) 

The  terminal  constraint  requires  that  the  missile  intercepts  the  target  at  the  final 
time.  The  state  constraint  is  given  by 

51  =a-20  +  /^i  =  0  (5.145) 

52  =  -«  +  /22  =  0  (5.146) 

which  implies  that  the  angle  of  attack  must  remain  between  0°  and  20°.  la  is  the 
slack  variable  added  to  the  a  state  since  the  constraint  is  an  inequality  constraint. 
The  control  constraints  are  expressed  as 

Cl  =«-20  +  /2j  =0  (5.147) 

C2  = -u-20  +  ll2  =  0.  (5.148) 

The  control  is  limited  to  a  value  of  -20  to  20  g’s.  For  simplification  define  the  state 
vector  as  x 

x  =  [x  y  h  V  7  V’  tr  a  WY  (5.149) 

then 

x  =  f{x,u,t)  (5.150) 

where  f(x,Ujt)  is  the  system  of  equations  given  above.  The  Bolza  function,  G,  and 
the  variational  Hamiltonian,  H,  are  expressed  as 

G{tf,Xf,v)  =  tf  +  v  (xf-  XT^y  +  (yf  -  yxjY  +  (hj  - 
H{t,x,u,\,$,  fi,l)  =  X^V  cos  7  cos  ^  -f-  XyV  cos  7  sin  ^  +  XhV  sin  7 
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+Av  (Tsina  —  jD)  —  ^sin7 

+A~  (T  sin  a  +  Z)  cos  cr  —  ^  cos  7 
IVW  V  . 


+A^  -777;^ - (T  sin  a  +  X)  sin  (T 

VW  cos  7 


+AcCc  {glrady  +  ^  ^  ■  (T sina  +  L) 

H — -  l^arctan  ^  _  crj  _ 

-\-9\  (a  —  20  +  +  02  (—a  +  /^2) 

+/I1  —  20  +  +  fi2  M  —  20  —  1^2'} 


(5.151) 


In  addition  to  satisfying  the  constraints,  the  remaining  necessary  conditions  are 
expressed  as 


Hf  =  -a 


'2v  (xf  -  XTf) 

2v  (yf  -  VTf) 

2v  (hf  ~  hxj) 


(5.152) 


Gl  = 


(5.153) 


Ax  —  0 


Aj,  —  0 


A.  =0 


Xy  =  —Xx  COS  7  cos  ij:  ~  Xy  cos  7  sin  ^  —  A/,  sin  7 


(5.154) 

(5.155) 

(5.156) 
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+A^  — — - (r  sin  a  +  L)  sin  <7 

[  K '“Vv  COS  7 

=  X^V  sin  7  cos  +  XyV  sin7sin^  —  XhV  cos  7 

q  ,  r  flsin7  . 

+Av5'cos7  -  A^— sin7  +  A^  - 5~(^ +  L)sma 

y  \  VW  cos^  7 

A^  =  A^y  cos  j  sin  if)  —  XyV  cos  7  cos  if) 

Ac  =  Ay:^(rsina  -  D)  -  X^^^{T  cos  a  +  ^)  cos  cr 

-X^p-r^ - (Tcosa  +  L)  sincT  +  XaCa^^^^-{T  cosa  +  L) 

VW  cos  7 

—0l  +  02 

O!  +  L)  sin  <7  -  sin  a  +  L)  cos  <7  + 

Xw  =  Av^(Tcosa-D)  +  A^:^^(Tsina  +  Z,)cost7 


(5.157) 


(5.158) 

(5.159) 


(T  sin  a  +  L)  cos  cr 
.  ^  phase. ^ 


+^V'T77i7| - sin  O'  +  X)  sin  <7  -  XcCc,-r^{T  sin  a  +  L) 

^  VW^  cos  7 


(5.160) 

^5.161) 

'TfT 


(5.162) 


Hi  = 


2/ii/ui 

2p2lu2  . 


/f.  =  0 


(5.163) 


(5.164) 


Examining  the  optimal  control  law,  i/„  =  0 


-  KCo,  {gLi.  +  ““)  ’  “  (  y  ;  „ 


+,..  +  ..  =  0 


(5.165) 


and  the  necessary  condition  Hi  =  0 


Hi  = 


2pi  lui 
2p2iu2  . 


(5.166) 
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three  cases  are  required  to  be  examined  depending  upon  the  value  of  u  to  satisfy 
Eqs  (5.166)  and  (5.165). 

CASE  1:  u  =  =  20.  This  case  occurs  when  the  constraint  at  u^ax  is 

active.  Hence,  is  zero  because  the  constraint  at  Umin  is  inactive.  /„i  is  zero  to 
satisfy  Eq.  (5.147)  and  1^2  =  to  satisfy  Eq.  (5.148).  An  expression  for  ni  is 
obtained  from  Eq.  (5.165)  with  the  substitution  u  =  20  yielding 


-20A^a  (; 


(5.167) 


CASE  2:  u  =  Umin  =  —  20.  This  case  occurs  when  the  constraint  at  Umin  is 
active.  Hence,  fii  is  zero  because  the  constraint  at  u^ax  is  inactive.  ly,2  is  zero  to 
satisfy  Eq.  (5.148)  and  =  ^/40  to  satisfy  Eq.  (5.147).  An  expression  for  fi2  is 
obtained  from  Eq.  (5.165)  with  the  substitution  u  —  —20  yielding 


fl2  —  20\aCo:  {glmdy  +  20^^  ^ 


Qcmdy 


Kdlmdy  +  202  ^ 


(5.168) 


CASE  3:  Umin  <  u  <  Umax-  This  case  occurs  when  u  is  between  the  upper  and 
lower  constraint  values.  Hence,  =  M2  =  0  because  both  constraints  are  inactive. 
lu\  =  V2O  —  u  to  satisfy  Eq.  (5.147)  and  =  \/20  +  u  to  satisfy  Eq.  (5.148). 
Equation  (5.165)  is  used  to  determine  the  value  of  u.  Adding  a  common  denominator, 
Eq.  (5.165)  becomes 

(glmdy  +  W  =  —9cmdy-  (5.169) 

'^<7 

Squaring  both  sides  of  the  equation  and  expanding  yields 

+  glmdy'^^  ~  glmdy  =  0 

where 

K  =  -  . 

T^XqCq, 


(5.171) 


(5.170) 
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Solving  for  u  produces  four  possible  solutions: 


'  {l  [-9cmdy  +  ^/gLdy+^^<^9cmdy]y 
-  {i  [~9cmdy  +  yjgLdy  +  ^K^gLdylY 

{I  {-gLdy  -  Ygtmdy+'^K^gLdylY 
■  -  {I  [~9Ldy  -  YgLdy  +  ^K^gLdylY 


The  last  two  possibilities  are  always  imaginary  and  can  immediately  be  discounted. 
Hence,  the  remaining  solutions  are 


{i  [-gLdy  +  YgLdy  +  ^K^9Ldy]Y 

-  {I  [-glndy  +  Y9tmdy-^^K‘^9lmdy\Y 


An  interesting  problem  arises  when  gcmdy  is  zero  which  implies  flight  is  occur¬ 
ring  in  a  vertical  plane.  Solving  for  u  with  gcmdy  =  0  indicates  u  =  0.  When  these 
values  are  substituted  back  into  Eq.  (5.165),  the  result  can  only  be  evaluated  in  the 
limit  as  u  goes  to  zero.  Additionally,  the  numerical  value  of  acmd  is  either  0°  or  180°, 
depending  upon  the  manner  in  which  u  approaches  zero.  If  u  approaches  zero  from 
below,  then  acmd  is  180°  and  if  u  approaches  zero  from  above,  then  acmd  is  Oo.  Thus, 
acmd  is  not  an  explicit  function  of  u.  Therefore,  Eq.  (5.164)  with  a  constant  value 
applied  for  acmd  reduces  to 


XaCa  =  0.  (5.174) 

An  expression  for  the  control  is  not  contained  in  Eq.  (5.174).  For  this  third  case, 
flight  in  the  vertical  plane  is  a  special  case  and  will  be  examined  in  more  detail. 

5.5  Two-Dimensional  Pursuit 

The  differential  equations  for  the  two-dimensional  pursuit  reduce  to 

X  =  V  cos  7  (5.175) 
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A  =  y  sin7 

(5.176) 

V  =  {T  cos  a  —  D)  —  gsin-y 

(5.177) 

7  =  ^^(rsma-f-T)S  ^cos7 

(5.178) 

a  =  Ca  u  —  —  {T  sin  a  -f  T)  S 

(5.179) 

1 

II 

(5.180) 

f  l,u  >  0 

E=  { 

(  — 1,M  <  0 

(5.181) 

where 


S  is  representing  the  roll  in  this  model.  Recall  the  roll  command  for  the  three- 
dimensional  intercept 

(^cmd  =  arctan  •  (5.182) 

When  gcmdy  is  zero  and  Ccmd  is  0°  for  positive  u  and  180°  for  negative  n  if  a  four 
quadrant  arctan  is  used.  The  terminal  constraints  reduce  to 


=  (®/  -  ^Tf)  +  {hf  -  hTf)  J  =  0. 


(5.183) 


The  state  and  control  constraints  remain  the  same  as  in  the  three-dimensional  pursuit 
and  are  given  by  Eqs  (5.145)  through  (5.148).  The  Bolza  and  Hamiltonian  functions 
are 

G  =  tf  +  v  (xf-XTj)  +{hf-hTj)^=0  (5.184) 

H  =  \xV  cos  7  -f  \hV  sin  7  -f  Ay  {T  cos  a  —  Z?)  —  sin  7 

+\aCa  u- ■^{Ts'ma  + L)'E^ 

—XwWf  -f  01  —  20  /^i^  -f-  02  (“Q:  -H  /^2) 

-f/^i  —  20  -|-  /„i^  -|-  H2  —  20  -(- .  (5.185) 
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The  necessary  conditions  are 


Aa;  —  0 


(5.186) 


A;.  =  0 

2g 

=  -As  cos  7  -  A/j  sin  7  - 

V^W  ^  ^ 

+A«C«— E 

=  X^V  sin  7  —  XhV  cos  7  +  Xyg  cos  7  —  X^-^  sin  7 
=  AK^Tsinc-  A,  (^Tcosa)  S 
+Ac,  cos  a^'L  + 01  + $2 

Aiv  =  Ay-j^  (T  cos  a  -  D)  +  a  +  X)  S 


-A,-j;j^(rsina  +  X)E 


=  -l 


’2y  (a:/  — 

2v  (hf  —  hxf^ 


(5.187) 

(5.188) 


(5.189) 

(5.190) 

(5.191) 

(5.192) 


Xf  =  Gi.  = 


(5.193) 


2^2^02 

Hi=  =0 

2/^1  lui 

_2^2^u2. 

Hu  —  XofCct  "f“  g\  /^2  0 

Hun  =  0. 


(5.194) 


(5.195) 

(5.196) 
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As  in  the  three-dimensional  pursuit,  three  cases  must  be  examined  to  satisfy  Eqs  (5.194) 
and  (5.195). 

CASE  1:  u  =  Umax  =  20.  This  case  occurs  when  the  constraint  at  Umax  is 
active.  Hence,  /i2  is  zero  because  the  constraint  at  Umin  is  inactive.  /„i  is  zero  to 
satisfy  Eq.  (5.147)  and  lu2  —  '\/40  to  satisfy  Eq.  (5.148).  An  expression  for  fii  is 
obtained  from  Eq.  (5.195)  yielding 


Ml  =  -KCa.  (5.197) 

CASE  2.  u  =  Umin  =  —20.  This  case  occurs  when  the  constraint  at  Umin  is 
active.  Hence,  /ii  is  zero  because  the  constraint  at  Umax  is  inactive.  Iu2  is  zero  to 
satisfy  Eq.  (5.148)  and  lui  =  i/iO  to  satisfy  Eq.  (5.147).  An  expression  for  //2  is 
obtained  from  Eq.  (5.195)  yielding 


M2  =  ^aCa-  (5.198) 

CASE  3:  Umin  <  u  <  Umax-  This  case  occurs  when  u  is  between  the  upper  and 
lower  constraint  values.  Hence,  //^  =  /12  =  0  because  both  constraints  are  inactive. 
lui  =  \/20  —  u  to  satisfy  Eq.  (5.147)  and  1^2  =  ^20  +  u  to  satisfy  Eq.  (5.148). 
Examining  Eq.  (5.195)  to  determine  the  value  of  u  yields 

Aa  =  0.  (5.199) 

The  necessary  condition  is  satisfied  when  is  zero.  An  expression  for  u  is  not 
contained  in  the  optimality  condition.  However,  since  Huu  is  singular,  the  extremal 
arc  for  A^  =  0  is  called  a  singular  arc  (2:246).  It  is  still  possible  to  determine  an 
expression  for  u  for  some  finite  time  interval  that  satisfies  the  optimality  condition. 
The  time  rate  of  change  of  the  optimality  condition  must  remain  zero  along  the 
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singular  arc  (2:252).  Hence, 


I  (H^)  =  LCa  (5.200) 

=  Ca  Ay  sin  a  —  cos  S  +  A^C^  cos  S(  5.201) 

=  0.  (5.202) 

This  expression  still  does  not  contain  the  control,  so  the  second  time  derivative  of 

Hu  must  be  zero.  The  differential  equations  for  the  system  can  be  expressed  as 

x  =  f{x)+g{x)u  (5.203) 

where 

V  cos  7 

V  sin  7 

4r  (T  cos  a  —  D)  —  o  sin  7 

f=  ^  '  (5.204) 

(T  sin  q:  +  T)  S  -  cos  7 

-^(rsina  +  L)E 

—wj 

0  1 
0 
0 

(5.205) 

0 

nf 

0  . 

Note  that  Eq.  (5.200)  can  be  expressed  as 

~  {Hu)  =  {gj  -  f^g) 

=  \^q  (5.206) 
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where 


Q  =  9xf-  fx9 
0 

0 

Ca-^T  sin  a 

cos  q;S 
Ciscos  aS 

0 


With  the  above  definitions  for  f,  g,  and  q,  the  second  time  derivative  of 

1  can  be 

expressed  as 

(P  \  xT- 

—  (Hu)  =  X^q  +  X  q 

(5.208) 

=  ><^9x  if  +  9^)  -  X^  {f^  +  g^u)  q 

(5.209) 

=  i9xf  -  fx9)  +  {9x9  -  9x9)  u 

(5.210) 

=  0 

(5.211) 

Noting  that 

^0.  =0 

(5.212) 

f  =  x-gu 

(5.213) 

and  recalling  that  is  identically  zero  for  the  singular  arc,  an  expression  for  u  is 


{qj  -  f^q) 

(9x9) 

XxA  +  XhB  —  XyC  —  A-yZ) 
XyE  +  X^F 


(5.214) 

(5.215) 
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where 


A  =  Ca-^T  (sin  a  cos  7  +  cos  a  sin  7S) 
B  —  Ca-^T  (sin  a  sin  7  —  cos  a  cos  7S) 


c  =  c.ir 


-777  (i  sin  a  +  L)  cos  q;L  +  — 7  sin  a  +  —77  sin  a 
W  W  VW 


2Dg  . 


cos  aS  +  Coi—  sin  a  cos  aS 
V  W 


D  =  Ca 


VW  S--rcosasm7S 


T 

— -^T  (T  sin  a  +  L)  sin  a  —  cos  aE  +  y^T  (T  sin  a  +  L)  sinaS 

—  -~^LT  sinaS  —  sin  a  cos  7  +  cos  a  cos  7S  —  777 cos^a 
VW  V  V  W  . 

E  =  Cl-^T  cos  a 

F^Cl-^TsmaT. 

VW 


Having  determined  u  for  all  three  cases,  the  two-dimensional  pursuit  is  now  formu¬ 
lated  in  such  a  matter  that  the  shooting  method  previously  discussed  can  be  used 
to  find  the  optimal  trajectory. 
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VI.  Conclusions  and  Recommendation 


The  simulation  was  successfully  recreated  in  MATLAB^^.  Errors  discovered 
in  ENGAGE  were  analyzed  to  determine  the  impact  of  the  errors.  The  recreated 
simulation  will  permit  future  numerical  analysis  of  the  missile’s  performance.  The 
necessary  and  sufficient  conditions  for  an  extremal  solution  to  the  free  final  time 
optimal  control  problem  were  presented.  The  shooting  method  was  developed  to 
provide  a  numerical  solution  to  the  optimal  control  problem  posed  as  a  two-point 
boundary-value  problem.  An  example  problem  of  a  lunar  lift  off  was  formulated  as 
a  two-point  boundary-value  problem  and  then  numerically  solved  via  the  shooting 
method.  The  necessary  conditions  for  an  extremal  solution  were  applied  to  the  three- 
dimensional  pursuit.  The  vertical  plane  intercept  scenario  was  determined  to  be  an 
interesting  case.  Therefore,  the  two-dimensional  intercept  was  analyzed  and  a  two- 
point  boundary-value  problem  was  formulated  to  determine  the  optimal  trajectory 
for  a  minimum  time  target  intercept. 

Future  efforts  in  this  project  should  be  directed  towards  numerically  solving 
the  two-dimensional  pursuit  to  analyze  the  behavior  of  the  optimal  loft  command 
for  various  initial  conditions.  Efforts  should  then  expand  the  problem  to  the  three- 
dimensional  pursuit  to  develop  an  expression  for  the  loft  control  which  approximates 
the  optimal  loft  control. 
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